SCM

SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
ViewVC logotype

Diff of /pkg/Matrix/R/diagMatrix.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

pkg/R/diagMatrix.R revision 1331, Sat Jul 22 17:59:53 2006 UTC pkg/Matrix/R/diagMatrix.R revision 3079, Tue Mar 31 15:29:43 2015 UTC
# Line 1  Line 1 
1    #### All methods for "diagonalMatrix" and its subclasses,
2    ####  currently "ddiMatrix", "ldiMatrix"
3    
4  ## Purpose: Constructor of diagonal matrices -- ~= diag() ,  ## Purpose: Constructor of diagonal matrices -- ~= diag() ,
5  ##          but *not* diag() extractor!  ##          but *not* diag() extractor!
6  Diagonal <- function(n, x = NULL)  Diagonal <- function(n, x = NULL)
7  {  {
8      ## Allow  Diagonal(4)  and  Diagonal(x=1:5)      ## Allow  Diagonal(4), Diagonal(x=1:5), and  Diagonal(4, TRUE)
9      if(missing(n))      n <- if(missing(n)) length(x) else {
         n <- length(x)  
     else {  
10          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)          stopifnot(length(n) == 1, n == as.integer(n), n >= 0)
11          n <- as.integer(n)          as.integer(n)
12      }      }
13    
14      if(missing(x)) # unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
15          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
16      else {      else {
17          stopifnot(length(x) == n)          lx <- length(x)
18            lx.1 <- lx == 1L
19            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
20          if(is.logical(x))          if(is.logical(x))
21              cl <- "ldiMatrix"              cl <- "ldiMatrix"
22          else {          else if(is.numeric(x)) {
23              cl <- "ddiMatrix"              cl <- "ddiMatrix"
24              x <- as.numeric(x)              x <- as.numeric(x)
25          }          }
26          new(cl, Dim = c(n,n), diag = "N", x = x)          else if(is.complex(x)) {
27                cl <- "zdiMatrix"  # will not yet work
28            } else stop("'x' has invalid data type")
29            if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
30                new(cl, Dim = c(n,n), diag = "U")
31            else
32                new(cl, Dim = c(n,n), diag = "N",
33                    x = if(lx.1) rep.int(x,n) else x)
34      }      }
35  }  }
36    
37  setAs("diagonalMatrix", "triangularMatrix",  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38        function(from) {                              shape = if(missing(cols)) "t" else "g",
39            n <- from@Dim[1]                              unitri, kind,
40            i <- seq(length = n)                              cols = if(n) 0:(n - 1L) else integer(0))
41            x <- from@x  {
42            new(if(is.numeric(x)) "dtTMatrix" else "ltTMatrix",      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
43        if(!(mcols <- missing(cols)))
44            stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
45        m <- length(cols)
46        if(missing(kind))
47            kind <-
48                if(is.double(x)) "d"
49                else if(is.logical(x)) "l"
50                else { ## for now
51                    storage.mode(x) <- "double"
52                    "d"
53                }
54        else stopifnot(any(kind == c("d","l","n")))
55        stopifnot(is.character(shape), nchar(shape) == 1,
56                  any(shape == c("t","s","g"))) # triangular / symmetric / general
57        if((missing(unitri) || unitri) && shape == "t" &&
58           (mcols || cols == 0:(n-1L)) &&
59           ((any(kind == c("l", "n")) && allTrue(x)) ||
60            (    kind == "d"          && allTrue(x == 1)))) { ## uni-triangular
61            new(paste0(kind,"tCMatrix"), Dim = c(n,n),
62                       uplo = uplo, diag = "U", p = rep.int(0L, n+1L))
63        }
64        else if(kind == "n") {
65            if(shape == "g")
66                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
67            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
68                     i = cols, p = 0:m)
69        }
70        else { ## kind != "n" -- have x slot :
71            if((lx <- length(x)) == 1) x <- rep.int(x, m)
72            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
73            if(shape == "g")
74                new(paste0(kind, "gCMatrix"), Dim = c(n,m),
75                    x = x, i = cols, p = 0:m)
76            else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
77                     x = x, i = cols, p = 0:m)
78        }
79    }
80    
81    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83        .sparseDiagonal(n, x, uplo, shape = "s")
84    
85    # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86    .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87        .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
88    
89    
90    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
91    ## Bert's code built on a post by Andy Liaw who most probably was influenced
92    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
93    ## who posted his bdiag() function written in December 1995.
94    if(FALSE)##--- no longer used:
95    .bdiag <- function(lst) {
96        ## block-diagonal matrix [a dgTMatrix] from list of matrices
97        stopifnot(is.list(lst), length(lst) >= 1)
98        dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
99        ## make sure we had all matrices:
100        if(!(is.matrix(dims) && nrow(dims) == 2))
101            stop("some arguments are not matrices")
102        csdim <- rbind(rep.int(0L, 2),
103                       apply(dims, 1, cumsum))
104        r <- new("dgTMatrix")
105        r@Dim <- as.integer(csdim[nrow(csdim),])
106        add1 <- matrix(1:0, 2,2)
107        for(i in seq_along(lst)) {
108            indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
109            if(is.null(dim(indx))) ## non-square matrix
110                r[indx[[1]],indx[[2]]] <- lst[[i]]
111            else ## square matrix
112                r[indx[,1], indx[,2]] <- lst[[i]]
113        }
114        r
115    }
116    ## expand(<mer>) needed something like bdiag() for lower-triangular
117    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
118    ##  implementation for those; now extended and generalized:
119    .bdiag <- function(lst) {
120        ## block-diagonal matrix [a dgTMatrix] from list of matrices
121        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123        Tlst <- lapply(lapply(lst, as_Csp2), # includes "diagU2N"
124                       as, "TsparseMatrix")
125        if(nl == 1) return(Tlst[[1]])
126        ## else
127        i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
128        j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
129    
130        clss <- vapply(Tlst, class, "")
131        typ <- substr(clss, 2, 2)
132        knd <- substr(clss, 1, 1)
133        sym <- typ == "s" # symmetric ones
134        tri <- typ == "t" # triangular ones
135        use.n <- any(is.n <- knd == "n")
136        if(use.n && !(use.n <- all(is.n))) {
137            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
138            knd [is.n] <- "l"
139        }
140        use.l <- !use.n && all(knd == "l")
141        if(all(sym)) { ## result should be *symmetric*
142            uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
143            tLU <- table(uplos)# of length 1 or 2 ..
144            if(length(tLU) == 1) { ## all "U" or all "L"
145                useU <- uplos[1] == "U"
146            } else { ## length(tLU) == 2, counting "L" and "U"
147                useU <- diff(tLU) >= 0
148                if(useU && (hasL <- tLU[1] > 0))
149                    Tlst[hasL] <- lapply(Tlst[hasL], t)
150                else if(!useU && (hasU <- tLU[2] > 0))
151                    Tlst[hasU] <- lapply(Tlst[hasU], t)
152            }
153            if(use.n) { ## return nsparseMatrix :
154                r <- new("nsTMatrix")
155            } else {
156                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
157                r@x <- unlist(lapply(Tlst, slot, "x"))
158            }
159            r@uplo <- if(useU) "U" else "L"
160        }
161        else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
162                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
163           ){ ## *triangular* result
164    
165            if(use.n) { ## return nsparseMatrix :
166                r <- new("ntTMatrix")
167            } else {
168                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
169                r@x <- unlist(lapply(Tlst, slot, "x"))
170            }
171            r@uplo <- ULs[1L]
172        }
173        else {
174            if(any(sym))
175                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
176            if(use.n) { ## return nsparseMatrix :
177                r <- new("ngTMatrix")
178            } else {
179                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
180                r@x <- unlist(lapply(Tlst, slot, "x"))
181            }
182        }
183        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
184        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
185        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
186        r
187    }
188    
189    bdiag <- function(...) {
190        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
191        if(nA == 1 && !is.list(...))
192            return(as(..., "CsparseMatrix"))
193        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
194        if(length(alis) == 1)
195            return(as(alis[[1]], "CsparseMatrix"))
196    
197        ## else : two or more arguments
198        as(.bdiag(alis), "CsparseMatrix")
199    }
200    
201    setMethod("tril", "diagonalMatrix", function(x, k = 0, ...)
202        if(k >= 0) x else .setZero(x, paste0(.M.kind(x), "tCMatrix")))
203    
204    setMethod("triu", "diagonalMatrix", function(x, k = 0, ...)
205        if(k <= 0) x else  .setZero(x, paste0(.M.kind(x), "tCMatrix")))
206    
207    
208    
209    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
210        ## to triangular Tsparse
211        i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
212        new(paste0(kind, "tTMatrix"),
213                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,                diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
214                x = x, i = i, j = i)          uplo = uplo,
215            })          x = from@x, # <- ok for diag = "U" and "N" (!)
216            i = i, j = i)
217    }
218    
219  setAs("diagonalMatrix", "matrix",  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
220        function(from) {      ## to symmetric Tsparse
221            n <- from@Dim[1]            n <- from@Dim[1]
222            diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE      i <- seq_len(n) - 1L
223                                       } else from@x,      new(paste0(kind, "sTMatrix"),
224                 nrow = n, ncol = n)          Dim = from@Dim, Dimnames = from@Dimnames,
225        })          i = i, j = i, uplo = uplo,
226            x = if(from@diag == "N") from@x else ## "U"-diag
227            rep.int(switch(kind,
228                           "d" = 1.,
229                           "l" =,
230                           "n" = TRUE,
231                           ## otherwise
232                           stop(gettextf("%s kind not yet implemented",
233                                         sQuote(kind)), domain=NA)),
234                    n))
235    }
236    
237  setAs("diagonalMatrix", "generalMatrix",  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
238    diag2tT.u <- function(d, x, kind = .M.kind(d))
239        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
240    
241    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
242    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
243        clx <- getClassDef(class(x))
244        if(extends(clx, "symmetricMatrix"))
245            .diag2sT(d, uplo = x@uplo, kind)
246        else
247            .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
248    }
249    
250    ## FIXME: should not be needed {when ddi* is dsparse* etc}:
251    setMethod("is.finite", signature(x = "diagonalMatrix"),
252              function(x) is.finite(.diag2tT(x)))
253    setMethod("is.infinite", signature(x = "diagonalMatrix"),
254              function(x) is.infinite(.diag2tT(x)))
255    
256    ## In order to evade method dispatch ambiguity warnings,
257    ## and because we can save a .M.kind() call, we use this explicit
258    ## "hack"  instead of signature  x = "diagonalMatrix" :
259    ##
260    ## ddi*:
261    di2tT <- function(from) .diag2tT(from, "U", "d")
262    setAs("ddiMatrix", "triangularMatrix", di2tT)
263    ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
264    ## needed too (otherwise <dense> -> Tsparse is taken):
265    setAs("ddiMatrix", "TsparseMatrix", di2tT)
266    setAs("ddiMatrix", "dsparseMatrix", di2tT)
267    setAs("ddiMatrix", "CsparseMatrix",
268          function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
269    setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
270    ##
271    ## ldi*:
272    di2tT <- function(from) .diag2tT(from, "U", "l")
273    setAs("ldiMatrix", "triangularMatrix", di2tT)
274    ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
275    ## needed too (otherwise <dense> -> Tsparse is taken):
276    setAs("ldiMatrix", "TsparseMatrix", di2tT)
277    setAs("ldiMatrix", "lsparseMatrix", di2tT)
278    setAs("ldiMatrix", "CsparseMatrix",
279          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280    setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281    rm(di2tT)
282    
283    setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
285            x <- as(from, "matrix")            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            as(x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287               if(is.logical(x)) "lgeMatrix"                Dim = from@Dim, Dimnames = from@Dimnames)
 ## Not yet:  
 ##              else if(is.complex(x)) "zgeMatrix"  
 ##              else if(is.integer(x)) "igeMatrix"  
              else "dgeMatrix")  
288        })        })
289    
290  setAs("ddiMatrix", "dgTMatrix",  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
       function(from) {  
           n <- from@Dim[1]  
           i <- seq(length = n) - 1:1  
           new("dgTMatrix", i = i, j = i,  
               x = if(from@diag == "U") rep(1,n) else from@x,  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
291    
292  setAs("ddiMatrix", "dgCMatrix",  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293        function(from) as(as(from, "dgTMatrix"), "dgCMatrix"))  mkDiag <- function(x, n) {
294        y <- matrix(as0(mod=mode(x)), n,n)
295        if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
296        y
297    }
298    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
299    ## --> using diag() where possible, ==> .ddi2mat()
300    
301  setAs("ldiMatrix", "lgTMatrix",  .diag2mat <- function(from)
302        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
303            n <- from@Dim[1]      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
304            i <- (if(from@diag == "U") seq(length = n) else which(from@x)) - 1:1  
305            new("lgTMatrix", i = i, j = i,  .ddi2mat <- function(from)
306                Dim = c(n,n), Dimnames = from@Dimnames) })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
307    
308  setAs("ldiMatrix", "lgCMatrix",  setAs("ddiMatrix", "matrix", .ddi2mat)
309        function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))  ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
310    setAs("diagonalMatrix", "matrix", .diag2mat)
311  setAs("diagonalMatrix", "sparseMatrix",  
312        function(from)  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
313            as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))            function(x, mode) {
314                  n <- x@Dim[1]
315                  mod.x <- mode(x@x)
316                  r <- vector(mod.x, length = n^2)
317                  if(n)
318                      r[1 + 0:(n - 1L) * (n + 1)] <-
319                          if(x@diag == "U") as1(mod=mod.x) else x@x
320                  r
321              })
322    
323    setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
324          function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
325    
326    setAs("diagonalMatrix", "denseMatrix",
327          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
328    
329    ..diag.x <- function(m)                   rep.int(as1(m@x), m@Dim[1])
330    .diag.x  <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
331    
332    .diag.2N <- function(m) {
333        if(m@diag == "U") m@diag <- "N"
334        m
335    }
336    
337  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
338        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
339    setAs("ddiMatrix", "ddenseMatrix",
340          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
341    setAs("ldiMatrix", "ldenseMatrix",
342          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
343    
344    
345  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
346        function(from) {        function(from) {
347            d <- dim(from)            d <- dim(from)
348            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
349            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
350                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to \"diagonalMatrix\"")
351            x <- diag(from)            x <- diag(from)
352            if(is.logical(x)) {            if(is.logical(x)) {
353                cl <- "ldiMatrix"                cl <- "ldiMatrix"
354                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
355            } else {            } else {
356                cl <- "ddiMatrix"                cl <- "ddiMatrix"
357                uni <- all(x == 1)                uni <- allTrue(x == 1)
358                storage.mode(x) <- "double"                storage.mode(x) <- "double"
359            }            } ## TODO: complex
360            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
361                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
362        })        })
# Line 110  Line 371 
371            x <- diag(from)            x <- diag(from)
372            if(is.logical(x)) {            if(is.logical(x)) {
373                cl <- "ldiMatrix"                cl <- "ldiMatrix"
374                uni <- all(x)                uni <- allTrue(x)
375            } else {            } else {
376                cl <- "ddiMatrix"                cl <- "ddiMatrix"
377                uni <- all(x == 1)                uni <- allTrue(x == 1)
378                storage.mode(x) <- "double"                storage.mode(x) <- "double"
379            }            } ## TODO: complex
380            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
381                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
382        })        })
383    
 setMethod("t", signature(x = "diagonalMatrix"),  
           function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })  
   
 setMethod("isDiagonal", signature(object = "diagonalMatrix"),  
           function(object) TRUE)  
 setMethod("isTriangular", signature(object = "diagonalMatrix"),  
           function(object) TRUE)  
 setMethod("isSymmetric", signature(object = "diagonalMatrix"),  
           function(object) TRUE)  
384    
385  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
386            function(x = 1, nrow, ncol = n) {            function(x = 1, nrow, ncol) .diag.x(x))
387               if(x@diag == "U")  
388                   rep.int(if(is.logical(x@x)) TRUE else 1, x@Dim[1])  subDiag <- function(x, i, j, ..., drop) {
389               else x@x      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
390        x <- if(missing(i))
391            x[, j, drop=drop]
392        else if(missing(j))
393            if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
394        else
395            x[i,j, drop=drop]
396        if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
397    }
398    
399    setMethod("[", signature(x = "diagonalMatrix", i = "index",
400                             j = "index", drop = "logical"), subDiag)
401    setMethod("[", signature(x = "diagonalMatrix", i = "index",
402                             j = "missing", drop = "logical"),
403              function(x, i, j, ..., drop) {
404                  na <- nargs()
405                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
406                  if(na == 4)
407                       subDiag(x, i=i, , drop=drop)
408                  else subDiag(x, i=i,   drop=drop)
409              })
410    setMethod("[", signature(x = "diagonalMatrix", i = "missing",
411                             j = "index", drop = "logical"),
412              function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
413    
414    ## When you assign to a diagonalMatrix, the result should be
415    ## diagonal or sparse ---
416    ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
417    ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
418    replDiag <- function(x, i, j, ..., value) {
419        x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
420        if(missing(i))
421            x[, j] <- value
422        else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
423            na <- nargs()
424    ##         message("diagnosing replDiag() -- nargs()= ", na)
425            if(na == 4)
426                x[i, ] <- value
427            else if(na == 3)
428                x[i] <- value
429            else stop(gettextf("Internal bug: nargs()=%d; please report",
430                               na), domain=NA)
431        } else
432            x[i,j] <- value
433        if(isDiagonal(x)) as(x, "diagonalMatrix") else x
434    }
435    
436    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
437                                    j = "index", value = "replValue"), replDiag)
438    
439    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
440                                    j = "missing", value = "replValue"),
441                     function(x,i,j, ..., value) {
442                         ## message("before replDiag() -- nargs()= ", nargs())
443                         if(nargs() == 3)
444                             replDiag(x, i=i, value=value)
445                         else ## nargs() == 4 :
446                             replDiag(x, i=i, , value=value)
447            })            })
448    
449  setMethod("!", "ldiMatrix", function(e1) {  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
450      if(e1@diag == "N")                                  j = "index", value = "replValue"),
451          e1@x <- !e1@x                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
452      else { ## "U"  
453          e1@diag <- "N"  ## x[] <- value :
454          e1@x <- rep.int(FALSE, e1@Dim[1])  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
455                                    j = "missing", value = "ANY"),
456                     function(x,i,j, ..., value)
457                 {
458                  if(all0(value)) { # be faster
459                      r <- new(paste0(.M.kindC(getClassDef(class(x))),"tTMatrix"))# of all "0"
460                      r@Dim <- x@Dim
461                      r@Dimnames <- x@Dimnames
462                      r
463                  } else { ## typically non-sense: assigning to full sparseMatrix
464                      x[TRUE] <- value
465                      x
466      }      }
467              })
468    
469    
470    setReplaceMethod("[", signature(x = "diagonalMatrix",
471                                    i = "matrix", # 2-col.matrix
472                                    j = "missing", value = "replValue"),
473                     function(x,i,j, ..., value) {
474                         if(ncol(i) == 2) {
475                             if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
476                                 if(x@diag == "U") {
477                                     one <- as1(x@x)
478                                     if(any(value != one | is.na(value))) {
479                                         x@diag <- "N"
480                                         x@x <- rep.int(one, x@Dim[1])
481                                     } else return(x)
482                                 }
483                                 x@x[ii] <- value
484                                 x
485                             } else { ## no longer diagonal, but remain sparse:
486                                 x <- as(x, "TsparseMatrix")
487                                 x[i] <- value
488      x      x
489                             }
490                         }
491                         else { # behave as "base R": use as if vector
492                             x <- as(x, "matrix")
493                             x[i] <- value
494                             Matrix(x)
495                         }
496  })  })
497    
 ## Basic Matrix Multiplication {many more to add}  
498    
499  ## FIXME: extend this for 'ldi', i.e. do "diagonalMatrix"  ## value = "sparseMatrix":
500    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
501                                    value = "sparseMatrix"),
502                     function (x, i, j, ..., value)
503                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
504    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
505                                    value = "sparseMatrix"),
506                     function (x, i, j, ..., value)
507                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
508    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
509                                    value = "sparseMatrix"),
510                     function (x, i, j, ..., value)
511                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
512    
513    ## value = "sparseVector":
514    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
515                                    value = "sparseVector"),
516                     replDiag)
517    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
518                                    value = "sparseVector"),
519                     replDiag)
520    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
521                                    value = "sparseVector"),
522                     replDiag)
523    
524    
525    setMethod("t", signature(x = "diagonalMatrix"),
526              function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
527    
528    setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
529    setMethod("isTriangular", "diagonalMatrix", function(object, upper=NA, ...) TRUE)
530    setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
531    
532    setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
533    setMethod("skewpart", signature(x = "diagonalMatrix"), function(x) .setZero(x))
534    
535    setMethod("chol", signature(x = "ddiMatrix"),
536              function(x, pivot, ...) {
537                  if(x@diag == "U") return(x)
538                  ## else
539                  if(any(x@x < 0))
540                      stop("chol() is undefined for diagonal matrix with negative entries")
541                  x@x <- sqrt(x@x)
542                  x
543              })
544    ## chol(L) is L for logical diagonal:
545    setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
546    
547    setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
548              function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
549    
550    setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
551              function(x, type, ...) {
552                  if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
553                  type <- toupper(substr(type[1], 1, 1))
554                  isU <- (x@diag == "U") # unit-diagonal
555                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
556                  else { ## norm == "I","1","O","M" :
557                      if(isU) 1 else max(abs(x@x))
558                  }
559              })
560    
561    
562    
563    ## Basic Matrix Multiplication {many more to add}
564    ##       ---------------------
565    ## Note that "ldi" logical are treated as numeric
566  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
567      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      dimCheck(x,y)
568      if(x@diag != "U") {      if(x@diag != "U") {
569          if(y@diag != "U") x@x <- x@x * y@x          if(y@diag != "U") {
570                nx <- x@x * y@x
571                if(is.numeric(nx) && !is.numeric(x@x))
572                    x <- as(x, "dMatrix")
573                x@x <- as.numeric(nx)
574            }
575          return(x)          return(x)
576      } else ## x is unit diagonal      } else ## x is unit diagonal
577      return(y)      return(y)
578  }  }
579    
580  setMethod("%*%", signature(x = "ddiMatrix", y = "ddiMatrix"),  ##' Boolean Algebra/Arithmetic Product of Diagonal Matrices
581            diagdiagprod, valueClass = "ddiMatrix")  ##'  %&%
582    diagdiagprodBool <- function(x, y) {
583        dimCheck(x,y)
584        if(x@diag != "U") {
585            if(!is.logical(x@x)) x <- as(x, "lMatrix")
586            if(y@diag != "U") {
587                nx <- x@x & y@x
588                x@x <- as.logical(nx)
589            }
590            x
591        } else { ## x is unit diagonal: return y
592            if(!is.logical(y@x)) y <- as(y, "lMatrix")
593            y
594        }
595    }
596    
597  formals(diagdiagprod) <- alist(x=, y=NULL)  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
 setMethod("crossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  
           diagdiagprod, valueClass = "ddiMatrix")  
 setMethod("tcrossprod", signature(x = "ddiMatrix", y = "ddiMatrix"),  
598            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
599    
600    setMethod("%&%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
601              diagdiagprodBool, valueClass = "ldiMatrix")# do *not* have "ndiMatrix" !
602    
603    ##' Both Numeric or Boolean Algebra/Arithmetic Product of Diagonal Matrices
604    diagdiagprodFlexi <- function(x, y=NULL, boolArith = NA, ...)
605    {
606        dimCheck(x,y)
607        bool <- isTRUE(boolArith)
608        if(x@diag != "U") {
609            if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
610            if(y@diag != "U") {
611                if(bool) {
612                    nx <- x@x & y@x
613                    x@x <- as.logical(nx)
614                } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
615                    nx <- x@x * y@x
616                    if(is.numeric(nx) && !is.numeric(x@x))
617                        x <- as(x, "dMatrix")
618                    x@x <- as.numeric(nx)
619                }
620            }
621            x
622        } else { ## x is unit diagonal: return y
623            if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
624            y
625        }
626    }
627    setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
628              diagdiagprodFlexi)
629    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
630              diagdiagprodFlexi)
631    
632    ##' crossprod(x) := x'x
633    diagprod <- function(x, y = NULL, boolArith = NA, ...) {
634        bool <- isTRUE(boolArith)
635        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
636        if(x@diag != "U") {
637            if(bool) {
638                nx <- x@x & y@x
639                x@x <- as.logical(nx)
640            } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
641                nx <- x@x * x@x
642                if(is.numeric(nx) && !is.numeric(x@x))
643                    x <- as(x, "dMatrix")
644                x@x <- as.numeric(nx)
645            }
646        }
647        x
648    }
649    setMethod( "crossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
650    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
651    
652    
653    ## analogous to matdiagprod() below:
654  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
655      dx <- dim(x)      ## x is diagonalMatrix
656      dy <- dim(y)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
657      if(dx[2] != dy[1]) stop("non-matching dimensions")      Matrix(if(x@diag == "U") y else x@x * y)
     n <- dx[1]  
     as(if(x@diag == "U") y else x@x * y, "Matrix")  
658  }  }
659    setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
660    
661  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  formals(diagmatprod) <- alist(x=, y=NULL, boolArith = NA, ...=) ## FIXME boolArith
662            diagmatprod)  diagmatprod2 <- function(x, y=NULL, boolArith = NA, ...) {
663  formals(diagmatprod) <- alist(x=, y=NULL)      ## x is diagonalMatrix
664  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
665            diagmatprod)      Matrix(if(x@diag == "U") y else x@x * y)
666    }
667    setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod2)
668    
669  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
670      dx <- dim(x)      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
     dy <- dim(y)  
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
671      if(x@diag != "U")      if(x@diag != "U")
672          y@x <- x@x * y@x          y@x <- x@x * y@x
673      y      y
674  }  }
675  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
676            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
677  formals(diagdgeprod) <- alist(x=, y=NULL)  
678  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  diagGeprodBool <- function(x, y) {
679            diagdgeprod, valueClass = "dgeMatrix")      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
680        if(!is.logical(y@x)) y <- as(y, "lMatrix")
681        if(x@diag != "U")
682            y@x <- x@x & y@x
683        y
684    }
685    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
686    
687    diagGeprod2 <- function(x, y=NULL, boolArith = NA, ...) {
688        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
689        bool <- isTRUE(boolArith)
690        if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
691        if(x@diag != "U")
692            y@x <- if(bool) x@x & y@x else x@x * y@x
693        y
694    }
695    setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), diagGeprod2)
696    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), diagGeprod2)
697    
698  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  
699            function(x, y) {  ## analogous to diagmatprod() above:
700    matdiagprod <- function(x, y) {
701                dx <- dim(x)                dx <- dim(x)
702                dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
703                if(dx[2] != dy[1]) stop("non-matching dimensions")      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
704                as(if(y@diag == "U") x else x * rep.int(y@x, dx[1]), "Matrix")  }
705            })  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
706    
707  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
708                dx <- dim(x)                dx <- dim(x)
709                dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
               if(dx[2] != dy[1]) stop("non-matching dimensions")  
710                if(y@diag == "N")                if(y@diag == "N")
711                    x@x <- x@x * rep.int(y@x, dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
712                x                x
713    }
714    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
715    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
716    
717    gediagprodBool <- function(x, y) {
718        dx <- dim(x)
719        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
720        if(!is.logical(x@x)) x <- as(x, "lMatrix")
721        if(y@diag == "N")
722            x@x <- x@x & rep(y@x, each = dx[1])
723        x
724    }
725    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
726    
727    setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"),
728              function(x, y=NULL, boolArith = NA, ...) {
729                  dx <- dim(x)
730                  if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
731                  bool <- isTRUE(boolArith)
732                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
733                  Matrix(if(y@diag == "U") x else
734                         if(bool) x & rep(y@x, each = dx[1])
735                         else     x * rep(y@x, each = dx[1]))
736            })            })
737    
738    setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
739              function(x, y=NULL, boolArith = NA, ...) {
740                  dx <- dim(x)
741                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
742                  bool <- isTRUE(boolArith)
743                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
744                  Matrix(if(bool) t(rep.int(y@x, dx[2]) & x)
745                             else t(rep.int(y@x, dx[2]) * x))
746              })
747    
748    
749    gediagprod2 <- function(x, y=NULL, boolArith = NA, ...) {
750        dx <- dim(x)
751        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
752        bool <- isTRUE(boolArith)
753        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
754        if(y@diag == "N")
755            x@x <- if(bool) x@x & rep(y@x, each = dx[1])
756                   else     x@x * rep(y@x, each = dx[1])
757        x
758    }
759    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), gediagprod2)
760    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), gediagprod2)
761    
762    
763  ## crossprod {more of these}  ## crossprod {more of these}
764    
765  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
766    
767    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
768              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
769    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
770              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
771    
772    
773  ## FIXME:  ## FIXME:
774  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
775  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
776  ##           })  ##           })
777    
778  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  ##' @param x CsparseMatrix
779  ##        function(x, y = NULL) {  ##' @param y diagonalMatrix
780  ##           })  ##' @return x %*% y
781    Cspdiagprod <- function(x, y, boolArith = NA, ...) {
782        if((m <- ncol(x)) != y@Dim[1]) stop("non-matching dimensions")
783        if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
784            x <- .Call(Csparse_diagU2N, x)
785            cx <- getClass(class(x))
786            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
787                x <- as(x, "generalMatrix")
788            ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
789            if(isTRUE(boolArith)) {
790                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
791                x@x <- r <- x@x & y@x[x@i + 1L]
792                if(!anyNA(r) && !extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
793            } else {
794                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
795                x@x <- x@x * y@x[ind]
796            }
797            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
798                ## instead of dropping all factors, be smart about some
799                ## TODO ......
800                x@factors <- list()
801            }
802            x
803        } else { #  y is unit-diagonal ==> "return x"
804            cx <- getClass(class(x))
805            if(isTRUE(boolArith)) {
806                is.l <- if(extends(cx, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
807                    x <- as(x, "lMatrix"); TRUE } else extends(cx, "lMatrix")
808                if(is.l && !anyNA(x@x)) as(drop0(x), "nMatrix")
809                else if(is.l) x else # defensive:
810                as(x, "lMatrix")
811            } else {
812                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
813                if(extends(cx, "nMatrix") || extends(cx, "lMatrix"))
814                    as(x, "dMatrix") else x
815            }
816        }
817    }
818    
819    ##' @param x diagonalMatrix
820    ##' @param y CsparseMatrix
821    ##' @return x %*% y
822    diagCspprod <- function(x, y, boolArith = NA, ...) {
823        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
824        if(x@diag == "N") {
825            y <- .Call(Csparse_diagU2N, y)
826            cy <- getClass(class(y))
827            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
828                y <- as(y, "generalMatrix")
829            if(isTRUE(boolArith)) {
830                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
831                y@x <- r <- y@x & x@x[y@i + 1L]
832                if(!anyNA(r) && !extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
833            } else {
834                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
835                y@x <- y@x * x@x[y@i + 1L]
836            }
837            if(.hasSlot(y, "factors") && length(y@factors)) {
838         ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
839                ## instead of dropping all factors, be smart about some
840                ## keep <- character()
841                ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
842                ##     keep <- "LU"
843                ## }
844                ## y@factors <- yf[keep]
845                y@factors <- list()
846            }
847            y
848        } else { ## x @ diag  == "U"
849            cy <- getClass(class(y))
850            if(isTRUE(boolArith)) {
851                is.l <- if(extends(cy, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
852                    y <- as(y, "lMatrix"); TRUE } else extends(cy, "lMatrix")
853                if(is.l && !anyNA(y@x)) as(drop0(y), "nMatrix")
854                else if(is.l) y else # defensive:
855                as(y, "lMatrix")
856            } else {
857                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
858                if(extends(cy, "nMatrix") || extends(cy, "lMatrix"))
859                    as(y, "dMatrix") else y
860            }
861        }
862    }
863    
864  ### ---------------- diagonal  o   sparse  -----------------------------  ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
865    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
866              function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, y, boolArith=boolArith))
867    
868  ## These are cheap implementations via coercion  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
869              function(x, y=NULL, boolArith=NA, ...)
870                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
871    
872    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
873    ##  x'y == (y'x)'
874    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
875              function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, x, boolArith=boolArith)))
876    
877    setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
878              function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
879    
880    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
881              function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(y), boolArith=boolArith))
882    
883    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
884              function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
885    
886    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
887              function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(x, y, boolArith=boolArith))
888    
889    setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
890              function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
891    
892    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
893              function(x, y) diagCspprod(x, y, boolArith=NA))
894    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
895              function(x, y) diagCspprod(x, y, boolArith=TRUE))
896    
897  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
898    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
899    
900  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
901            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
902    
903  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
904            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
905    
906  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
907            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
908    
909  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
910            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
911    
912  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  }
           function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })  
913    
914  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
915            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y) Cspdiagprod(x, y, boolArith=NA))
916    setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
917              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
918    
919    ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
920    ##       do indeed work by going through sparse (and *not* ddense)!
921    
922    setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
923              function(a, b, ...) {
924                  a@x <- 1/ a@x
925                  a@Dimnames <- a@Dimnames[2:1]
926                  a
927              })
928    
929    solveDiag <- function(a, b, ...) {
930        if(a@Dim[1] != nrow(b))
931            stop("incompatible matrix dimensions")
932        ## trivially invert a 'in place' and multiply:
933        a@x <- 1/ a@x
934        a@Dimnames <- a@Dimnames[2:1]
935        a %*% b
936    }
937    setMethod("solve", signature(a = "diagonalMatrix", b = "matrix"),
938              solveDiag)
939    setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
940              solveDiag)
941    
942    ## Schur()  ---> ./eigen.R
943    
944    
945    
946    ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
947    
948    ## Use function for several signatures, in order to evade
949    diagOdiag <- function(e1,e2) {
950        ## result should also be diagonal _ if possible _
951        r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
952        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
953        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
954                           if(is.numeric(e2@x)) 0 else FALSE)
955        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
956            if(is.numeric(r)) { # "double" *or* "integer"
957                if(is.numeric(e2@x)) {
958                    e2@x <- r; return(.diag.2N(e2)) }
959                if(!is.numeric(e1@x))
960                    ## e.g. e1, e2 are logical;
961                    e1 <- as(e1, "dMatrix")
962                if(!is.double(r)) r <- as.double(r)
963            }
964            else if(is.logical(r))
965                e1 <- as(e1, "lMatrix")
966            else stop(gettextf("intermediate 'r' is of type %s",
967                               typeof(r)), domain=NA)
968            e1@x <- r
969            .diag.2N(e1)
970        }
971        else { ## result not diagonal, but at least symmetric:
972            ## e.g., m == m
973            isNum <- (is.numeric(r) || is.numeric(r00))
974            isLog <- (is.logical(r) || is.logical(r00))
975            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
976            d <- e1@Dim
977            n <- d[1]
978            stopifnot(length(r) == n)
979            if(isNum && !is.double(r)) r <- as.double(r)
980            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
981            newcl <-
982                paste0(if(isNum) "d" else if(isLog) {
983                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
984                } else stop("not yet implemented .. please report"), "syMatrix")
985    
986            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
987        }
988    }
989    
990    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
991    ## we use this hack instead of signature  x = "diagonalMatrix" :
992    diCls <- names(getClass("diagonalMatrix")@subclasses)
993    if(FALSE) {
994    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
995              diagOdiag)
996    } else { ## These are just for method disambiguation:
997        for(c1 in diCls)
998            for(c2 in diCls)
999                setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
1000    }
1001    
1002    ## diagonal  o  triangular  |-->  triangular
1003    ## diagonal  o  symmetric   |-->  symmetric
1004    ##    {also when other is sparse: do these "here" --
1005    ##     before conversion to sparse, since that loses "diagonality"}
1006    diagOtri <- function(e1,e2) {
1007        ## result must be triangular
1008        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
1009        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
1010        e1.0 <- if(is.numeric(d1)) 0 else FALSE
1011        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
1012        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
1013            diag(e2) <- r
1014            ## check what happens "in the triangle"
1015            e2.2 <- if(.n2) 2 else TRUE
1016            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
1017                n <- dim(e2)[1L]
1018                it <- indTri(n, upper = (e2@uplo == "U"))
1019                e2[it] <- callGeneric(e1.0, e2[it])
1020            }
1021            e2
1022        }
1023        else { ## result not triangular ---> general
1024            rr <- as(e2, "generalMatrix")
1025            diag(rr) <- r
1026            rr
1027        }
1028    }
1029    
1030    
1031    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
1032              diagOtri)
1033    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
1034    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
1035    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1036              function(e1,e2)
1037          { ## this must only trigger for *dense* e1
1038              switch(.Generic,
1039                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
1040                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
1041                     "*" = {
1042                         n <- e2@Dim[1L]
1043                         d2 <- if(e2@diag == "U") { # unit-diagonal
1044                             d <- rep.int(as1(e2@x), n)
1045                             e2@x <- d
1046                             e2@diag <- "N"
1047                             d
1048                         } else e2@x
1049                         e2@x <- diag(e1) * d2
1050                         e2
1051                     },
1052                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
1053                         e1 ^ as(e2, "denseMatrix")
1054                     },
1055                     ## otherwise:
1056                     callGeneric(e1, diag2Tsmart(e2,e1)))
1057    })
1058    
1059    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
1060    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1061              .Cmp.swap)
1062    ## '&' and "|'  are commutative:
1063    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1064              function(e1,e2) callGeneric(e2, e1))
1065    
1066    ## For almost everything else, diag* shall be treated "as sparse" :
1067    ## These are cheap implementations via coercion
1068    
1069    ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
1070    ## and because we can save an .M.kind() call, we use this explicit
1071    ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
1072    ##
1073    ## ddi*:
1074    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
1075              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1076    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
1077              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1078    ## ldi*
1079    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
1080              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1081    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
1082              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1083    
1084    ## Ops:  Arith  --> numeric : "dMatrix"
1085    ##       Compare --> logical
1086    ##       Logic   --> logical: "lMatrix"
1087    
1088    ## Other = "numeric" : stay diagonal if possible
1089    ## ddi*: Arith: result numeric, potentially ddiMatrix
1090    for(arg2 in c("numeric","logical"))
1091    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
1092              function(e1,e2) {
1093                  n <- e1@Dim[1]
1094                  f0 <- callGeneric(0, e2)
1095                  if(all0(f0)) { # remain diagonal
1096                      L1 <- (le <- length(e2)) == 1L
1097                      if(e1@diag == "U") {
1098                          if(any((r <- callGeneric(1, e2)) != 1)) {
1099                              e1@diag <- "N"
1100                              e1@x[seq_len(n)] <- r # possibly recycling r
1101                          } ## else: result = e1  (is "U" diag)
1102                      } else {
1103                          r <- callGeneric(e1@x, e2)
1104                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1105                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1106                      }
1107                      e1
1108                  } else
1109                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
1110              })
1111    
1112    for(arg1 in c("numeric","logical"))
1113    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
1114              function(e1,e2) {
1115                  n <- e2@Dim[1]
1116                  f0 <- callGeneric(e1, 0)
1117                  if(all0(f0)) { # remain diagonal
1118                      L1 <- (le <- length(e1)) == 1L
1119                      if(e2@diag == "U") {
1120                          if(any((r <- callGeneric(e1, 1)) != 1)) {
1121                              e2@diag <- "N"
1122                              e2@x[seq_len(n)] <- r # possibly recycling r
1123                          } ## else: result = e2  (is "U" diag)
1124                      } else {
1125                          r <- callGeneric(e1, e2@x)
1126                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1127                          e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1128                      }
1129                      e2
1130                  } else
1131                      callGeneric(e1, diag2tT.u(e2,e1, "d"))
1132              })
1133    
1134    ## ldi* Arith --> result numeric, potentially ddiMatrix
1135    for(arg2 in c("numeric","logical"))
1136    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
1137              function(e1,e2) {
1138                  n <- e1@Dim[1]
1139                  f0 <- callGeneric(0, e2)
1140                  if(all0(f0)) { # remain diagonal
1141                      L1 <- (le <- length(e2)) == 1L
1142                      E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1143                      ## storage.mode(E@x) <- "double"
1144                      if(e1@diag == "U") {
1145                          if(any((r <- callGeneric(1, e2)) != 1)) {
1146                              E@diag <- "N"
1147                              E@x[seq_len(n)] <- r # possibly recycling r
1148                          } ## else: result = E  (is "U" diag)
1149                      } else {
1150                          r <- callGeneric(e1@x, e2)
1151                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1152                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1153                      }
1154                      E
1155                  } else
1156                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1157              })
1158    
1159    for(arg1 in c("numeric","logical"))
1160    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
1161              function(e1,e2) {
1162                  n <- e2@Dim[1]
1163                  f0 <- callGeneric(e1, 0)
1164                  if(all0(f0)) { # remain diagonal
1165                      L1 <- (le <- length(e1)) == 1L
1166                      E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1167                      ## storage.mode(E@x) <- "double"
1168                      if(e2@diag == "U") {
1169                          if(any((r <- callGeneric(e1, 1)) != 1)) {
1170                              E@diag <- "N"
1171                              E@x[seq_len(n)] <- r # possibly recycling r
1172                          } ## else: result = E  (is "U" diag)
1173                      } else {
1174                          r <- callGeneric(e1, e2@x)
1175                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1176                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1177                      }
1178                      E
1179                  } else
1180                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1181              })
1182    
1183    ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1184    ##
1185    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1186    if(FALSE) {
1187        selectMethod("<", c("numeric","lMatrix"))# Compare
1188        selectMethod("&", c("numeric","lMatrix"))# Logic
1189    } ## so we don't need to define a method here :
1190    
1191    for(arg2 in c("numeric","logical"))
1192    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1193              function(e1,e2) {
1194                  n <- e1@Dim[1]
1195                  f0 <- callGeneric(0, e2)
1196                  if(all0(f0)) { # remain diagonal
1197                      L1 <- (le <- length(e2)) == 1L
1198                      E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1199                      ## storage.mode(E@x) <- "logical"
1200                      if(e1@diag == "U") {
1201                          if(any((r <- callGeneric(1, e2)) != 1)) {
1202                              E@diag <- "N"
1203                              E@x[seq_len(n)] <- r # possibly recycling r
1204                          } ## else: result = E  (is "U" diag)
1205                      } else {
1206                          r <- callGeneric(e1@x, e2)
1207                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1208                          E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1209                      }
1210                      E
1211                  } else
1212                      callGeneric(diag2tT.u(e1,e2, "d"), e2)
1213              })
1214    
1215    ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1216    for(arg2 in c("numeric","logical"))
1217    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1218              function(e1,e2) {
1219                  n <- e1@Dim[1]
1220                  f0 <- callGeneric(FALSE, e2)
1221                  if(all0(f0)) { # remain diagonal
1222                      L1 <- (le <- length(e2)) == 1L
1223    
1224                      if(e1@diag == "U") {
1225                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1226                              e1@diag <- "N"
1227                              e1@x[seq_len(n)] <- r # possibly recycling r
1228                          } ## else: result = e1  (is "U" diag)
1229                      } else {
1230                          r <- callGeneric(e1@x, e2)
1231                          ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1232                          e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1233                      }
1234                      e1
1235                  } else
1236                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1237              })
1238    
1239    
1240    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1241    for(other in c("ANY", "Matrix", "dMatrix")) {
1242        ## ddi*:
1243        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1244                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1245        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1246                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1247        ## ldi*:
1248        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1249                  function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1250        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1251                  function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1252    }
1253    
1254    ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1255    if(FALSE) # now also contains "geMatrix"
1256    dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1257                           names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1258    dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1259    for(DI in diCls) {
1260        dMeth <- if(extends(DI, "dMatrix"))
1261            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1262        else # "lMatrix", the only other kind for now
1263            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1264        for(c2 in c(dense.subCl, "Matrix")) {
1265            for(Fun in c("*", "&")) {
1266                setMethod(Fun, signature(e1 = DI, e2 = c2),
1267                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1268                setMethod(Fun, signature(e1 = c2, e2 = DI),
1269                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1270            }
1271            setMethod("^", signature(e1 = c2, e2 = DI),
1272                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1273            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1274                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1275        }
1276    }
1277    
1278    ## Group methods "Math", "Math2" in                     --> ./Math.R
1279    
1280    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1281    ### ----------   the last 4: separately here
1282    for(cl in diCls) {
1283    setMethod("any", cl,
1284              function (x, ..., na.rm) {
1285                  if(any(x@Dim == 0)) FALSE
1286                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1287              })
1288    setMethod("all",  cl, function (x, ..., na.rm) {
1289        n <- x@Dim[1]
1290        if(n >= 2) FALSE
1291        else if(n == 0 || x@diag == "U") TRUE
1292        else all(x@x, ..., na.rm = na.rm)
1293    })
1294    setMethod("prod", cl, function (x, ..., na.rm) {
1295        n <- x@Dim[1]
1296        if(n >= 2) 0
1297        else if(n == 0 || x@diag == "U") 1
1298        else ## n == 1, diag = "N" :
1299            prod(x@x, ..., na.rm = na.rm)
1300    })
1301    
1302    setMethod("sum", cl,
1303              function(x, ..., na.rm) {
1304                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
1305                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
1306              })
1307    }
1308    
1309    ## The remaining ones are  max, min, range :
1310    
1311    setMethod("Summary", "ddiMatrix",
1312              function(x, ..., na.rm) {
1313                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
1314                  else if(x@diag == "U")
1315                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
1316                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
1317              })
1318    setMethod("Summary", "ldiMatrix",
1319              function(x, ..., na.rm) {
1320                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
1321                  else if(x@diag == "U")
1322                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
1323                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
1324              })
1325    
1326    
1327    
# Line 267  Line 1336 
1336      invisible(x)      invisible(x)
1337  }  }
1338    
1339    ## somewhat consistent with "print" for sparseMatrix :
1340    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1341    
1342  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1343            function(object) prDiag(object))            function(object) {
1344                  d <- dim(object)
1345                  cl <- class(object)
1346                  cat(sprintf('%d x %d diagonal matrix of class "%s"',
1347                              d[1], d[2], cl))
1348                  if(d[1] < 50) {
1349                      cat("\n")
1350                      prDiag(object)
1351                  } else {
1352                      cat(", with diagonal entries\n")
1353                      show(diag(object))
1354                      invisible(object)
1355                  }
1356              })
1357    
1358    rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1359       dense.subCl, diCls)# not used elsewhere
1360    
1361    setMethod("summary", signature(object = "diagonalMatrix"),
1362              function(object, ...) {
1363                  d <- dim(object)
1364                  r <- summary(object@x, ...)
1365                  attr(r, "header") <-
1366                      sprintf('%d x %d diagonal Matrix of class "%s"',
1367                              d[1], d[2], class(object))
1368                  ## use ole' S3 technology for such a simple case
1369                  class(r) <- c("diagSummary", class(r))
1370                  r
1371              })
1372    
1373    print.diagSummary <- function (x, ...) {
1374        cat(attr(x, "header"),"\n")
1375        class(x) <- class(x)[-1]
1376        print(x, ...)
1377        invisible(x)
1378    }

Legend:
Removed from v.1331  
changed lines
  Added in v.3079

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge