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

revision 2096, Fri Dec 7 17:44:44 2007 UTC revision 2207, Mon Jul 7 22:34:52 2008 UTC
# Line 16  Line 16 
16      if(missing(x)) ## unit diagonal matrix      if(missing(x)) ## unit diagonal matrix
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          stopifnot(length(x) == n)          lx <- length(x)
20            stopifnot(lx == 1 || lx == n) # but keep 'x' short for now
21          if(is.logical(x))          if(is.logical(x))
22              cl <- "ldiMatrix"              cl <- "ldiMatrix"
23          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 26  Line 27 
27          else if(is.complex(x)) {          else if(is.complex(x)) {
28              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
29          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
30          new(cl, Dim = c(n,n), diag = "N", x = x)          new(cl, Dim = c(n,n), diag = "N",
31                x = if(lx == 1) rep.int(x,n) else x)
32      }      }
33  }  }
34    
35    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
36    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {
37        stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
38        if((lx <- length(x)) == 1) x <- rep.int(x, n)
39        else if(lx != n) stop("length(x) must be 1 or n")
40        cls <-
41            if(is.double(x)) "dsCMatrix"
42            else if(is.logical(x)) "lsCMatrix"
43            else { ## for now
44                storage.mode(x) <- "double"
45                "dsCMatrix"
46            }
47        new(cls, Dim = c(n,n), x = x, uplo = uplo,
48            i = if(n) 0:(n - 1L) else integer(0), p = 0:n)
49    }
50    
51  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
52  ### Bert's code built on a post by Andy Liaw who most probably was influenced  ### Bert's code built on a post by Andy Liaw who most probably was influenced
53  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
54  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
55    if(FALSE)##--- currently unused:
56  bdiag <- function(...) {  .bdiag <- function(lst) {
57      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
58      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
59      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
     dims <- sapply(mlist, dim)  
60      ## make sure we had all matrices:      ## make sure we had all matrices:
61      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
62          stop("some arguments are not matrices")          stop("some arguments are not matrices")
63      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
64                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
65      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
66        r@Dim <- as.integer(csdim[nrow(csdim),])
67      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
68      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
69          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
70          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
71              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
72          else ## square matrix          else ## square matrix
73              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
74        }
75        r
76      }      }
77      ## slightly debatable if we really should return Csparse.. :  ### Doug Bates needed something like bdiag() for lower-triangular
78      as(ret, "CsparseMatrix")  ### (Tsparse) Matrices and provided a much more efficient implementation:
79    .bdiag <- function(lst) {
80        ### block-diagonal matrix [a dgTMatrix] from list of matrices
81        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
82    
83        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
84                       as, "TsparseMatrix")
85    
86        if(nl == 1) return(Tlst[[1]])
87        ## else
88        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
89        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
90        new("dgTMatrix", Dim = c(i_off[nl+1], j_off[nl + 1]),
91            i = unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k])),
92            j = unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])),
93            x = unlist(lapply(Tlst, slot, "x")))
94  }  }
95    
96  diag2tT <- function(from) {  bdiag <- function(...) {
97        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
98        if(nA == 1 && !is.list(...))
99            return(as(..., "CsparseMatrix"))
100        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
101        if(length(alis) == 1)
102            return(as(alis[[1]], "CsparseMatrix"))
103    
104        ## else : two or more arguments
105        as(.bdiag(alis), "CsparseMatrix")
106    }
107    
108    
109    .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
110        ## to triangular Tsparse
111      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L
112      new(paste(.M.kind(from), "tTMatrix", sep=''),      new(paste(kind, "tTMatrix", sep=''),
113          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
114            uplo = uplo,
115          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
116          i = i, j = i)          i = i, j = i)
117  }  }
118    
119  diag2sT <- function(from) { # to symmetric Tsparse  .diag2sT <- function(from, uplo = "U", kind = .M.kind(from)) {
120      i <- if(from@diag == "U") integer(0) else seq_len(from@Dim[1]) - 1L      ## to symmetric Tsparse
121      new(paste(.M.kind(from), "sTMatrix", sep=''),      n <- from@Dim[1]
122        i <- seq_len(n) - 1L
123        new(paste(kind, "sTMatrix", sep=''),
124          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
125          x = from@x, i = i, j = i)          i = i, j = i, uplo = uplo,
126            x = if(from@diag == "N") from@x else ## "U"-diag
127            rep.int(switch(kind,
128                           "d" = 1.,
129                           "l" =,
130                           "n" = TRUE,
131                           ## otherwise
132                           stop("'", kind,"' kind not yet implemented")), n))
133    }
134    
135    ## diagonal -> triangular,  upper / lower depending on "partner":
136    diag2tT.u <- function(d, x, kind = .M.kind(d))
137        .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
138    
139    ## diagonal -> sparse {triangular OR symmetric} (upper / lower) depending on "partner":
140    diag2Tsmart <- function(d, x, kind = .M.kind(d)) {
141        clx <- getClassDef(class(x))
142        if(extends(clx, "symmetricMatrix"))
143            .diag2sT(d, uplo = x@uplo, kind)
144        else
145            .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
146  }  }
147    
148  setAs("diagonalMatrix", "triangularMatrix", diag2tT)  
149  setAs("diagonalMatrix", "sparseMatrix", diag2tT)  ## In order to evade method dispatch ambiguity warnings,
150    ## and because we can save a .M.kind() call, we use this explicit
151    ## "hack"  instead of signature  x = "diagonalMatrix" :
152    ##
153    ## ddi*:
154    diag2tT <- function(from) .diag2tT(from, "U", "d")
155    setAs("ddiMatrix", "triangularMatrix", diag2tT)
156    setAs("ddiMatrix", "sparseMatrix", diag2tT)
157  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
158  setAs("diagonalMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
159  ## is better than this:  setAs("ddiMatrix", "CsparseMatrix",
160  ## setAs("diagonalMatrix", "sparseMatrix",        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
161  ##       function(from)  setAs("ddiMatrix", "symmetricMatrix",
162  ##        as(from, if(is(from, "dMatrix")) "dgCMatrix" else "lgCMatrix"))        function(from) .diag2sT(from, "U", "d"))
163  setAs("diagonalMatrix", "CsparseMatrix",  ##
164        function(from) as(diag2tT(from), "CsparseMatrix"))  ## ldi*:
165    diag2tT <- function(from) .diag2tT(from, "U", "l")
166    setAs("ldiMatrix", "triangularMatrix", diag2tT)
167    setAs("ldiMatrix", "sparseMatrix", diag2tT)
168    ## needed too (otherwise <dense> -> Tsparse is taken):
169    setAs("ldiMatrix", "TsparseMatrix", diag2tT)
170    setAs("ldiMatrix", "CsparseMatrix",
171          function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
172    setAs("ldiMatrix", "symmetricMatrix",
173          function(from) .diag2sT(from, "U", "l"))
174    
 setAs("diagonalMatrix", "symmetricMatrix", diag2sT)  
175    
176  setAs("diagonalMatrix", "matrix",  setAs("diagonalMatrix", "nMatrix",
177        function(from) {        function(from) {
178            n <- from@Dim[1]            n <- from@Dim[1]
179            diag(x = if(from@diag == "U") { if(is.numeric(from@x)) 1. else TRUE            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
180                                       } else from@x,            new("ntTMatrix", i = i, j = i, diag = from@diag,
181                 nrow = n, ncol = n)                Dim = from@Dim, Dimnames = from@Dimnames)
182        })        })
183    
 setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  
       function(from) as(from, paste(.M.kind(from), "gCMatrix", sep='')))  
   
 .diag.x <- function(m) {  
     if(m@diag == "U")  
         rep.int(if(is.numeric(m@x)) 1. else TRUE,  
                 m@Dim[1])  
     else m@x  
 }  
184    
185  .diag.2N <- function(m) {  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
186      if(m@diag == "U") m@diag <- "N"  mkDiag <- function(x, n) {
187      m      y <- matrix(as0(mod=mode(x)), n,n)
188        if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x
189        y
190  }  }
191    
192  ## given the above, the following  4  coercions should be all unneeded;  setAs("diagonalMatrix", "matrix",
 ## we prefer triangular to general:  
 setAs("ddiMatrix", "dgTMatrix",  
193        function(from) {        function(from) {
194            .Deprecated("as(, \"sparseMatrix\")")            ## want "ldiMatrix" -> <logical> "matrix" :
195            n <- from@Dim[1]            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
196            i <- seq_len(n) - 1L                   n = from@Dim[1])
197            new("dgTMatrix", i = i, j = i, x = .diag.x(from),        })
               Dim = c(n,n), Dimnames = from@Dimnames) })  
198    
199  setAs("ddiMatrix", "dgCMatrix",  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
200        function(from) as(as(from, "sparseMatrix"), "dgCMatrix"))            function(x, mode) {
201                  n <- x@Dim[1]
202                  mod.x <- mode(x@x)
203                  r <- vector(mod.x, length = n^2)
204                  if(n)
205                      r[1 + 0:(n - 1) * (n + 1)] <-
206                          if(x@diag == "U") as1(mod=mod.x) else x@x
207                  r
208              })
209    
210  setAs("ldiMatrix", "lgTMatrix",  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
211        function(from) {        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
           .Deprecated("as(, \"sparseMatrix\")")  
           n <- from@Dim[1]  
           if(from@diag == "U") { # unit-diagonal  
               x <- rep.int(TRUE, n)  
               i <- seq_len(n) - 1L  
           } else { # "normal"  
               nz <- nz.NA(from@x, na. = TRUE)  
               x <- from@x[nz]  
               i <- which(nz) - 1L  
           }  
           new("lgTMatrix", i = i, j = i, x = x,  
               Dim = c(n,n), Dimnames = from@Dimnames) })  
212    
213  setAs("ldiMatrix", "lgCMatrix",  .diag.x <- function(m) if(m@diag == "U") rep.int(as1(m@x), m@Dim[1]) else m@x
       function(from) as(as(from, "lgTMatrix"), "lgCMatrix"))  
214    
215    .diag.2N <- function(m) {
216        if(m@diag == "U") m@diag <- "N"
217        m
218    }
219    
 if(FALSE) # now have faster  "ddense" -> "dge"  
220  setAs("ddiMatrix", "dgeMatrix",  setAs("ddiMatrix", "dgeMatrix",
221        function(from) as(as(from, "matrix"), "dgeMatrix"))        function(from) .Call(dup_mMatrix_as_dgeMatrix, from))
222    setAs("ddiMatrix", "ddenseMatrix",
223          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
224    setAs("ldiMatrix", "ldenseMatrix",
225          function(from) as(as(from, "triangularMatrix"),"denseMatrix"))
226    
227    
228  setAs("matrix", "diagonalMatrix",  setAs("matrix", "diagonalMatrix",
229        function(from) {        function(from) {
# Line 188  Line 268 
268  setMethod("diag", signature(x = "diagonalMatrix"),  setMethod("diag", signature(x = "diagonalMatrix"),
269            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
270    
271    subDiag <- function(x, i, j, ..., drop) {
 subDiag <- function(x, i, j, drop) {  
272      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
273      x <- if(missing(i))      x <- if(missing(i))
274          x[, j, drop=drop]          x[, j, drop=drop]
# Line 197  Line 276 
276          x[i, , drop=drop]          x[i, , drop=drop]
277      else      else
278          x[i,j, drop=drop]          x[i,j, drop=drop]
279      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
280  }  }
281    
282  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
283                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
284  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
285                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
286            function(x, i, drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))
287  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
288                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
289            function(x, j, drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
290    
291  ## When you assign to a diagonalMatrix, the result should be  ## When you assign to a diagonalMatrix, the result should be
292  ## diagonal or sparse ---  ## diagonal or sparse ---
293  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
294  replDiag <- function(x, i, j, value) {  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
295    replDiag <- function(x, i, j, ..., value) {
296      x <- as(x, "sparseMatrix")      x <- as(x, "sparseMatrix")
     message(sprintf("diagnosing replDiag() -- nargs() = %d\n", nargs()))  
297      if(missing(i))      if(missing(i))
298          x[, j] <- value          x[, j] <- value
299      else if(missing(j))      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
300            na <- nargs()
301    ##         message("diagnosing replDiag() -- nargs()= ", na)
302            if(na == 4)
303          x[i, ] <- value          x[i, ] <- value
304      else          else if(na == 3)
305                x[i] <- value
306            else stop("Internal bug: nargs()=",na,"; please report")
307        } else
308          x[i,j] <- value          x[i,j] <- value
309      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
310  }  }
311    
312  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
313                                  j = "index", value = "replValue"), replDiag)                                  j = "index", value = "replValue"), replDiag)
314    
315  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index",
316                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
317                   function(x, i, value) replDiag(x, i=i, value=value))                   function(x,i,j, ..., value) {
318  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "matrix", # 2-col.matrix                       ## message("before replDiag() -- nargs()= ", nargs())
319                         if(nargs() == 3)
320                             replDiag(x, i=i, value=value)
321                         else ## nargs() == 4 :
322                             replDiag(x, i=i, , value=value)
323                     })
324    
325    setReplaceMethod("[", signature(x = "diagonalMatrix",
326                                    i = "matrix", # 2-col.matrix
327                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
328                   function(x, i, value) {                   function(x,i,j, ..., value) {
329                       if(ncol(i) == 2) {                       if(ncol(i) == 2) {
330                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only                           if(all((ii <- i[,1]) == i[,2])) { # replace in diagonal only
331                                 if(x@diag == "U") {
332                                     one <- as1(x@x)
333                                     if(any(value != one | is.na(value))) {
334                                         x@diag <- "N"
335                                         x@x <- rep.int(one, x@Dim[1])
336                                     }
337                                 }
338                               x@x[ii] <- value                               x@x[ii] <- value
339                               x                               x
340                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
# Line 251  Line 352 
352    
353  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
354                                  j = "index", value = "replValue"),                                  j = "index", value = "replValue"),
355                   function(x, j, value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
356    
357    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
358                                    value = "scarceMatrix"),
359                     function (x, i, j, ..., value)
360                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
361    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
362                                    value = "scarceMatrix"),
363                     function (x, i, j, ..., value)
364                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
365    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
366                                    value = "scarceMatrix"),
367                     function (x, i, j, ..., value)
368                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
369    
370    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
371                                    value = "sparseVector"),
372                     replDiag)
373    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
374                                    value = "sparseVector"),
375                     replDiag)
376    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
377                                    value = "sparseVector"),
378                     replDiag)
379    
380    
381  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
# Line 262  Line 386 
386  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isTriangular", signature(object = "diagonalMatrix"),
387            function(object) TRUE)            function(object) TRUE)
388  setMethod("isSymmetric", signature(object = "diagonalMatrix"),  setMethod("isSymmetric", signature(object = "diagonalMatrix"),
389            function(object) TRUE)            function(object, ...) TRUE)
390    
391  setMethod("chol", signature(x = "ddiMatrix"),# pivot = "ANY"  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
392            function(x, pivot) {  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)
393                if(any(x@x < 0)) stop("chol() is undefined for diagonal matrix with negative entries")  
394    setMethod("chol", signature(x = "ddiMatrix"),
395              function(x, pivot, ...) {
396                  if(x@diag == "U") return(x)
397                  ## else
398                  if(any(x@x < 0))
399                      stop("chol() is undefined for diagonal matrix with negative entries")
400                x@x <- sqrt(x@x)                x@x <- sqrt(x@x)
401                x                x
402            })            })
403  ## chol(L) is L for logical diagonal:  ## chol(L) is L for logical diagonal:
404  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot) x)  setMethod("chol", signature(x = "ldiMatrix"), function(x, pivot, ...) x)
405    
406    setMethod("determinant", signature(x = "diagonalMatrix", logarithm = "logical"),
407              function(x, logarithm, ...) mkDet(.diag.x(x), logarithm))
408    
409    setMethod("norm", signature(x = "diagonalMatrix", type = "character"),
410              function(x, type, ...) {
411                  if((n <- x@Dim[1]) == 0) return(0) # as for "sparseMatrix"
412                  type <- toupper(substr(type[1], 1, 1))
413                  isU <- (x@diag == "U") # unit-diagonal
414                  if(type == "F") sqrt(if(isU) n else sum(x@x^2))
415                  else { ## norm == "I","1","O","M" :
416                      if(isU) 1 else max(abs(x@x))
417                  }
418              })
419    
420    
421    
422  ## Basic Matrix Multiplication {many more to add}  ## Basic Matrix Multiplication {many more to add}
423  ##       ---------------------  ##       ---------------------
424  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
425  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
426      if(any(dim(x) != dim(y))) stop("non-matching dimensions")      n <- dimCheck(x,y)[1]
427      if(x@diag != "U") {      if(x@diag != "U") {
428          if(y@diag != "U") {          if(y@diag != "U") {
429              nx <- x@x * y@x              nx <- x@x * y@x
# Line 305  Line 451 
451    
452    
453  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
454        ## x is diagonalMatrix
455      dx <- dim(x)      dx <- dim(x)
456      dy <- dim(y)      dy <- dim(y)
457      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
458      n <- dx[1]      n <- dx[1]
459      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
460  }  }
   
461  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
462            diagmatprod)            diagmatprod)
463    ## sneaky .. :
464  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
465  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),
466            diagmatprod)            diagmatprod)
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
467    
468  diagdgeprod <- function(x, y) {  diagGeprod <- function(x, y) {
469      dx <- dim(x)      dx <- dim(x)
470      dy <- dim(y)      dy <- dim(y)
471      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
# Line 328  Line 473 
473          y@x <- x@x * y@x          y@x <- x@x * y@x
474      y      y
475  }  }
476  setMethod("%*%", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
477            diagdgeprod, valueClass = "dgeMatrix")  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
478  formals(diagdgeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL)
479  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
480            diagdgeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
481    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
482              diagGeprod)
483    
484  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  matdiagprod <- function(x, y) {
           function(x, y) {  
485                dx <- dim(x)                dx <- dim(x)
486                dy <- dim(y)                dy <- dim(y)
487                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
488                as(if(y@diag == "U") x else x * rep(y@x, each = dx[1]), "Matrix")      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
489            })  }
490    setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),
491              matdiagprod)
492    formals(matdiagprod) <- alist(x=, y=NULL)
493    setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),
494              matdiagprod)
495    
496  setMethod("%*%", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprod <- function(x, y) {
           function(x, y) {  
497                dx <- dim(x)                dx <- dim(x)
498                dy <- dim(y)                dy <- dim(y)
499                if(dx[2] != dy[1]) stop("non-matching dimensions")                if(dx[2] != dy[1]) stop("non-matching dimensions")
500                if(y@diag == "N")                if(y@diag == "N")
501                    x@x <- x@x * rep(y@x, each = dx[1])                    x@x <- x@x * rep(y@x, each = dx[1])
502                x                x
503            })  }
504    setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
505    setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
506    formals(gediagprod) <- alist(x=, y=NULL)
507    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
508              gediagprod)
509    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
510              gediagprod)
511    
512  ## crossprod {more of these}  ## crossprod {more of these}
513    
514  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
515    
516    setMethod("%*%", signature(x = "diagonalMatrix", y = "denseMatrix"),
517              function(x, y) if(x@diag == "U") y else x %*% as(y, "generalMatrix"))
518    setMethod("%*%", signature(x = "denseMatrix", y = "diagonalMatrix"),
519              function(x, y) if(y@diag == "U") x else as(x, "generalMatrix") %*% y)
520    
521    
522  ## FIXME:  ## FIXME:
523  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),  ## setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "denseMatrix"),
524  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
# Line 366  Line 529 
529  ##           })  ##           })
530    
531  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
532            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))
533    
534  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
535            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))
536    
537  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
538            function(x, y = NULL) { x <- as(x, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))
539    
540  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
541            function(x, y = NULL) { y <- as(y, "sparseMatrix"); callGeneric() })            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))
542    
543    
544  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()
545  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
546            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) as(x, "sparseMatrix") %*% y)
547    setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
548              function(x, y) x %*% as(y, "sparseMatrix"))
549  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)
550  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.
551  ## ==> do this:  ## ==> do this:
# Line 392  Line 557 
557  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
558  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
559    
 setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  
           function(x, y) x %*% as(y, "sparseMatrix"))  
560    
561    
562  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
# Line 416  Line 579 
579  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),  setMethod("solve", signature(a = "diagonalMatrix", b = "Matrix"),
580            solveDiag)            solveDiag)
581    
582    ## Schur()  ---> ./eigen.R
583    
584    
585    
586  ### ---------------- diagonal  o  sparse  -----------------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
   
587    
588  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
589  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)
590  diagOdiag <- function(e1,e2) { # result should also be diagonal  diagOdiag <- function(e1,e2) {
591        ## result should also be diagonal _ if possible _
592      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
593        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
594        r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
595                           if(is.numeric(e2@x)) 0 else FALSE)
596        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
597      if(is.numeric(r)) {      if(is.numeric(r)) {
598          if(is.numeric(e2@x)) {          if(is.numeric(e2@x)) {
599              e2@x <- r; return(.diag.2N(e2)) }              e2@x <- r; return(.diag.2N(e2)) }
# Line 439  Line 607 
607      e1@x <- r      e1@x <- r
608      .diag.2N(e1)      .diag.2N(e1)
609  }  }
610        else { ## result not diagonal, but at least symmetric:
611            isNum <- (is.numeric(r) || is.numeric(r00))
612            isLog <- (is.logical(r) || is.logical(r00))
613    
614            if(getOption("verbose"))
615                message("exploding  <diag>  o  <diag>  into dense matrix")
616            d <- e1@Dim
617            n <- d[1]
618            stopifnot(length(r) == n)
619            xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
620            newcl <-
621                paste(if(isNum) "d" else if(isLog) {
622                    if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
623                } else stop("not yet implemented .. please report")
624                      ,
625                      "syMatrix", sep='')
626    
627            new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
628        }
629    }
630    
631    ### This would be *the* way, but we get tons of "ambiguous method dispatch"
632    ## we use this hack instead of signature  x = "diagonalMatrix" :
633    diCls <- names(getClass("diagonalMatrix")@subclasses)
634    if(FALSE) {
635  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "diagonalMatrix"),
636            diagOdiag)            diagOdiag)
637  ## These two are just for method disambiguation:  } else { ## These are just for method disambiguation:
638  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "diagonalMatrix"),      for(c1 in diCls)
639            diagOdiag)          for(c2 in diCls)
640  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ddiMatrix"),              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
641            diagOdiag)  }
642    
643  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## FIXME:    diagonal  o  triangular  |-->  triangular
644  ## -----     diagonal  o  symmetric   |-->  symmetric  ## -----     diagonal  o  symmetric   |-->  symmetric
# Line 456  Line 648 
648  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
649  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
650    
651  ## for disambiguation  ## For disambiguation --- define this for "sparseMatrix" , then for "ANY";
652  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "sparseMatrix"),  ## and because we can save an .M.kind() call, we use this explicit
653            function(e1,e2) callGeneric(as(e1, "sparseMatrix"), e2))  ## "hack" for all diagonalMatrix *subclasses* instead of just "diagonalMatrix" :
654  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "diagonalMatrix"),  ##
655            function(e1,e2) callGeneric(e1, as(e2, "sparseMatrix")))  ## ddi*:
656  ## in general:  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "sparseMatrix"),
657  setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "ANY"),            function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
658            function(e1,e2) callGeneric(as(e1,"sparseMatrix"), e2))  setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ddiMatrix"),
659  setMethod("Ops", signature(e1 = "ANY", e2 = "diagonalMatrix"),            function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
660            function(e1,e2) callGeneric(e1, as(e2,"sparseMatrix")))  ## ldi*
661    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "sparseMatrix"),
662              function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
663    setMethod("Ops", signature(e1 = "sparseMatrix", e2 = "ldiMatrix"),
664              function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
665    
666    ## Ops:  Arith  --> numeric : "dMatrix"
667    ##       Compare --> logical
668    ##       Logic   --> logical: "lMatrix"
669    
670    ##  other = "numeric" : stay diagonal if possible
671    ## ddi*: Arith: result numeric, potentially ddiMatrix
672    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),
673              function(e1,e2) {
674                  n <- e1@Dim[1]; nsq <- n^2
675                  f0 <- callGeneric(0, e2)
676                  if(all(is0(f0))) { # remain diagonal
677                      L1 <- (le <- length(e2)) == 1L
678                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
679                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
680                          e1@diag <- "N"
681                          if(L1) r <- rep.int(r, n)
682                      } else
683                          r <- callGeneric(e1@x, e2)
684                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
685                      return(e1)
686                  }
687                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
688              })
689    
690    setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),
691              function(e1,e2) {
692                  n <- e2@Dim[1]; nsq <- n^2
693                  f0 <- callGeneric(e1, 0)
694                  if(all(is0(f0))) { # remain diagonal
695                      L1 <- (le <- length(e1)) == 1L
696                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
697                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
698                          e2@diag <- "N"
699                          if(L1) r <- rep.int(r, n)
700                      } else
701                          r <- callGeneric(e1, e2@x)
702                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
703                      return(e2)
704                  }
705                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
706              })
707    
708    ## ldi* Arith --> result numeric, potentially ddiMatrix
709    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),
710              function(e1,e2) {
711                  n <- e1@Dim[1]; nsq <- n^2
712                  f0 <- callGeneric(0, e2)
713                  if(all(is0(f0))) { # remain diagonal
714                      L1 <- (le <- length(e2)) == 1L
715                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
716                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
717                          e1@diag <- "N"
718                          if(L1) r <- rep.int(r, n)
719                      } else
720                          r <- callGeneric(e1@x, e2)
721                      e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))
722                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
723                      return(e1)
724                  }
725                  callGeneric(diag2tT.u(e1,e2, "d"), e2)
726              })
727    
728    setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),
729              function(e1,e2) {
730                  n <- e2@Dim[1]; nsq <- n^2
731                  f0 <- callGeneric(e1, 0)
732                  if(all(is0(f0))) { # remain diagonal
733                      L1 <- (le <- length(e1)) == 1L
734                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
735                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
736                          e2@diag <- "N"
737                          if(L1) r <- rep.int(r, n)
738                      } else
739                          r <- callGeneric(e1, e2@x)
740                      e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))
741                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
742                      return(e2)
743                  }
744                  callGeneric(e1, diag2tT.u(e2,e1, "d"))
745              })
746    
747    ## ddi*: for "Ops" without Arith --> result logical, potentially ldi
748    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),
749              function(e1,e2) {
750                  n <- e1@Dim[1]; nsq <- n^2
751                  f0 <- callGeneric(0, e2)
752                  if(all(is0(f0))) { # remain diagonal
753                      L1 <- (le <- length(e2)) == 1L
754                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
755                      if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {
756                          e1@diag <- "N"
757                          if(L1) r <- rep.int(r, n)
758                      } else
759                          r <- callGeneric(e1@x, e2)
760                      e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))
761                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
762                      return(e1)
763                  }
764                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
765              })
766    
767    setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),
768              function(e1,e2) {
769                  n <- e2@Dim[1]; nsq <- n^2
770                  f0 <- callGeneric(e1, 0)
771                  if(all(is0(f0))) { # remain diagonal
772                      L1 <- (le <- length(e1)) == 1L
773                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
774                      if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {
775                          e2@diag <- "N"
776                          if(L1) r <- rep.int(r, n)
777                      } else
778                          r <- callGeneric(e1, e2@x)
779                      e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))
780                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
781                      return(e2)
782                  }
783                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
784              })
785    
786    ## ldi*: for "Ops" without Arith --> result logical, potentially ldi
787    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),
788              function(e1,e2) {
789                  n <- e1@Dim[1]; nsq <- n^2
790                  f0 <- callGeneric(FALSE, e2)
791                  if(all(is0(f0))) { # remain diagonal
792                      L1 <- (le <- length(e2)) == 1L
793                      if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)
794                      if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {
795                          e1@diag <- "N"
796                          if(L1) r <- rep.int(r, n)
797                      } else
798                          r <- callGeneric(e1@x, e2)
799                      e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
800                      return(e1)
801                  }
802                  callGeneric(diag2tT.u(e1,e2, "l"), e2)
803              })
804    
805    setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),
806              function(e1,e2) {
807                  n <- e2@Dim[1]; nsq <- n^2
808                  f0 <- callGeneric(e1, FALSE)
809                  if(all(is0(f0))) { # remain diagonal
810                      L1 <- (le <- length(e1)) == 1L
811                      if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)
812                      if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {
813                          e2@diag <- "N"
814                          if(L1) r <- rep.int(r, n)
815                      } else
816                          r <- callGeneric(e1, e2@x)
817                      e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]
818                      return(e2)
819                  }
820                  callGeneric(e1, diag2tT.u(e2,e1, "l"))
821              })
822    
823    
824    
825    ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
826    for(other in c("ANY", "Matrix", "dMatrix")) {
827        ## ddi*:
828        setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
829                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))
830        setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
831                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))
832        ## ldi*:
833        setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
834                  function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))
835        setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
836                  function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))
837    }
838    ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as
839    ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:
840    for(cl in diCls) {
841        setMethod("&", signature(e1 = cl, e2 = "ANY"),
842                  function(e1,e2) e1 & as(e2,"Matrix"))
843        setMethod("&", signature(e1 = "ANY", e2 = cl),
844                  function(e1,e2) as(e1,"Matrix") & e2)
845        for(c2 in c("denseMatrix", "Matrix")) {
846            setMethod("&", signature(e1 = cl, e2 = c2),
847                      function(e1,e2) e1 & Diagonal(x = diag(e2)))
848            setMethod("&", signature(e1 = c2, e2 = cl),
849                      function(e1,e2) Diagonal(x = diag(e1)) & e2)
850        }
851    }
852    
853    
854    ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
855    ### ----------  any, all: separately here
856    for(cl in diCls) {
857    setMethod("any", cl,
858              function (x, ..., na.rm) {
859                  if(any(x@Dim == 0)) FALSE
860                  else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
861              })
862    setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))
863    setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))
864    
865    setMethod("sum", cl,
866              function(x, ..., na.rm) {
867                  r <- sum(x@x, ..., na.rm = na.rm)# double or integer, correctly
868                  if(x@diag == "U" && !is.na(r)) r + x@Dim[1] else r
869              })
870    }
871    
872    ## The remaining ones are  max, min, range :
873    
874    setMethod("Summary", "ddiMatrix",
875              function(x, ..., na.rm) {
876                  if(any(x@Dim == 0)) callGeneric(numeric(0), ..., na.rm=na.rm)
877                  else if(x@diag == "U")
878                      callGeneric(x@x, 0, 1, ..., na.rm=na.rm)
879                  else callGeneric(x@x, 0, ..., na.rm=na.rm)
880              })
881    setMethod("Summary", "ldiMatrix",
882              function(x, ..., na.rm) {
883                  if(any(x@Dim == 0)) callGeneric(logical(0), ..., na.rm=na.rm)
884                  else if(x@diag == "U")
885                      callGeneric(x@x, FALSE, TRUE, ..., na.rm=na.rm)
886                  else callGeneric(x@x, FALSE, ..., na.rm=na.rm)
887              })
888    
889    
890    
# Line 484  Line 903 
903            function(object) {            function(object) {
904                d <- dim(object)                d <- dim(object)
905                cl <- class(object)                cl <- class(object)
906                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
907                            d[1], d[2], cl))                            d[1], d[2], cl))
908                  if(d[1] < 50) {
909                      cat("\n")
910                prDiag(object)                prDiag(object)
911                  } else {
912                      cat(", with diagonal entries\n")
913                      show(diag(object))
914                      invisible(object)
915                  }
916            })            })

Legend:
Removed from v.2096  
changed lines
  Added in v.2207

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