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 2226, Mon Jul 21 17:15:17 2008 UTC pkg/Matrix/R/diagMatrix.R revision 3079, Tue Mar 31 15:29:43 2015 UTC
# Line 5  Line 5 
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          lx <- length(x)          lx <- length(x)
18          stopifnot(lx == 1 || lx == n) # but keep 'x' short for now          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 if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 27  Line 26 
26          else if(is.complex(x)) {          else if(is.complex(x)) {
27              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
28          } else stop("'x' has invalid data type")          } 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",          new(cl, Dim = c(n,n), diag = "N",
33              x = if(lx == 1) rep.int(x,n) else x)                  x = if(lx.1) rep.int(x,n) else x)
34      }      }
35  }  }
36    
37  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38                                shape = if(missing(cols)) "t" else "g",
39                                unitri, kind,
40                                cols = if(n) 0:(n - 1L) else integer(0))
41    {
42      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
43      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!(mcols <- missing(cols)))
44      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
45      stopifnot(is.character(shape), nchar(shape) == 1,      m <- length(cols)
46                any(shape == c("t","s","g"))) # triangular / symmetric / general      if(missing(kind))
47      kind <-      kind <-
48          if(is.double(x)) "d"          if(is.double(x)) "d"
49          else if(is.logical(x)) "l"          else if(is.logical(x)) "l"
# Line 45  Line 51 
51              storage.mode(x) <- "double"              storage.mode(x) <- "double"
52              "d"              "d"
53          }          }
54      new(paste(kind, shape, "CMatrix", sep=''),      else stopifnot(any(kind == c("d","l","n")))
55          Dim = c(n,n), x = x, uplo = uplo,      stopifnot(is.character(shape), nchar(shape) == 1,
56          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)                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()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83      .sparseDiagonal(n, x, uplo, shape = "s")      .sparseDiagonal(n, x, uplo, shape = "s")
84    
85  ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87      .sparseDiagonal(n, x, uplo, shape = "t")      .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.  ## 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  ## 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  ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
93  ### who posted his bdiag() function written in December 1995.  ## who posted his bdiag() function written in December 1995.
94  if(FALSE)##--- no longer used:  if(FALSE)##--- no longer used:
95  .bdiag <- function(lst) {  .bdiag <- function(lst) {
96      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
97      stopifnot(is.list(lst), length(lst) >= 1)      stopifnot(is.list(lst), length(lst) >= 1)
98      dims <- sapply(lst, dim, USE.NAMES=FALSE)      dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
99      ## make sure we had all matrices:      ## make sure we had all matrices:
100      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
101          stop("some arguments are not matrices")          stop("some arguments are not matrices")
# Line 92  Line 120 
120      ## block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
121      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"      Tlst <- lapply(lapply(lst, as_Csp2), # includes "diagU2N"
124                     as, "TsparseMatrix")                     as, "TsparseMatrix")
125      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
126      ## else      ## else
127      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))      i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
128      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))      j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
129    
130      clss <- sapply(Tlst, class)      clss <- vapply(Tlst, class, "")
131      knds <- substr(clss, 2, 2)      typ <- substr(clss, 2, 2)
132      sym  <- knds == "s" # symmetric ones      knd <- substr(clss, 1, 1)
133      tri  <- knds == "t" # triangular ones      sym <- typ == "s" # symmetric ones
134      use.n <- any(is.n <- substr(clss,1,1) == "n")      tri <- typ == "t" # triangular ones
135      if(use.n && !(use.n <- all(is.n)))      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")          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*      if(all(sym)) { ## result should be *symmetric*
142          uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"          uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
143          tLU <- table(uplos)# of length 1 or 2 ..          tLU <- table(uplos)# of length 1 or 2 ..
144          if(length(tLU) == 1) { ## all "U" or all "L"          if(length(tLU) == 1) { ## all "U" or all "L"
145              useU <- uplos[1] == "U"              useU <- uplos[1] == "U"
# Line 121  Line 153 
153          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
154              r <- new("nsTMatrix")              r <- new("nsTMatrix")
155          } else {          } else {
156              r <- new("dsTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
157              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
158          }          }
159          r@uplo <- if(useU) "U" else "L"          r@uplo <- if(useU) "U" else "L"
160      }      }
161      else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"      else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
162                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
163         ){ ## *triangular* result         ){ ## *triangular* result
164    
165          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
166              r <- new("ntTMatrix")              r <- new("ntTMatrix")
167          } else {          } else {
168              r <- new("dtTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
169              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
170          }          }
171          r@uplo <- ULs[1L]          r@uplo <- ULs[1L]
# Line 144  Line 176 
176          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
177              r <- new("ngTMatrix")              r <- new("ngTMatrix")
178          } else {          } else {
179              r <- new("dgTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
180              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
181          }          }
182      }      }
# Line 166  Line 198 
198      as(.bdiag(alis), "CsparseMatrix")      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)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
210      ## to triangular Tsparse      ## to triangular Tsparse
211      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
212      new(paste(kind, "tTMatrix", sep=''),      new(paste0(kind, "tTMatrix"),
213          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
214          uplo = uplo,          uplo = uplo,
215          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
# Line 181  Line 220 
220      ## to symmetric Tsparse      ## to symmetric Tsparse
221      n <- from@Dim[1]      n <- from@Dim[1]
222      i <- seq_len(n) - 1L      i <- seq_len(n) - 1L
223      new(paste(kind, "sTMatrix", sep=''),      new(paste0(kind, "sTMatrix"),
224          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
225          i = i, j = i, uplo = uplo,          i = i, j = i, uplo = uplo,
226          x = if(from@diag == "N") from@x else ## "U"-diag          x = if(from@diag == "N") from@x else ## "U"-diag
# Line 190  Line 229 
229                         "l" =,                         "l" =,
230                         "n" = TRUE,                         "n" = TRUE,
231                         ## otherwise                         ## otherwise
232                         stop("'", kind,"' kind not yet implemented")), n))                         stop(gettextf("%s kind not yet implemented",
233                                         sQuote(kind)), domain=NA)),
234                    n))
235  }  }
236    
237  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
238  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
239      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
240    
# Line 206  Line 247 
247          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)          .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,  ## In order to evade method dispatch ambiguity warnings,
257  ## and because we can save a .M.kind() call, we use this explicit  ## and because we can save a .M.kind() call, we use this explicit
258  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
259  ##  ##
260  ## ddi*:  ## ddi*:
261  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
262  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
263  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
264  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
265  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
266    setAs("ddiMatrix", "dsparseMatrix", di2tT)
267  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
268        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
269  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
       function(from) .diag2sT(from, "U", "d"))  
270  ##  ##
271  ## ldi*:  ## ldi*:
272  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
273  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
274  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
275  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
276  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
277    setAs("ldiMatrix", "lsparseMatrix", di2tT)
278  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
279        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
282    
283  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
           n <- from@Dim[1]  
285            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            new("ntTMatrix", i = i, j = i, diag = from@diag,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
288        })        })
289    
290    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
291    
292  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
294      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
295      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
296      y      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("diagonalMatrix", "matrix",  .diag2mat <- function(from)
302        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
303            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
304            mkDiag(if(from@diag == "U") as1(from@x) else from@x,  
305                   n = from@Dim[1])  .ddi2mat <- function(from)
306        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
307    
308    setAs("ddiMatrix", "matrix", .ddi2mat)
309    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
310    setAs("diagonalMatrix", "matrix", .diag2mat)
311    
312  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
313            function(x, mode) {            function(x, mode) {
# Line 263  Line 315 
315                mod.x <- mode(x@x)                mod.x <- mode(x@x)
316                r <- vector(mod.x, length = n^2)                r <- vector(mod.x, length = n^2)
317                if(n)                if(n)
318                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1L) * (n + 1)] <-
319                        if(x@diag == "U") as1(mod=mod.x) else x@x                        if(x@diag == "U") as1(mod=mod.x) else x@x
320                r                r
321            })            })
# Line 271  Line 323 
323  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
324        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        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  .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) {  .diag.2N <- function(m) {
# Line 291  Line 347 
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            } ## 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",
# Line 315  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        })        })
# Line 330  Line 386 
386            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
387    
388  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
389      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
390      x <- if(missing(i))      x <- if(missing(i))
391          x[, j, drop=drop]          x[, j, drop=drop]
392      else if(missing(j))      else if(missing(j))
393          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
394      else      else
395          x[i,j, drop=drop]          x[i,j, drop=drop]
396      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 344  Line 400 
400                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
401  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
402                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
403            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            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",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
411                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
412            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))            function(x, i, j, ..., drop) subDiag(x, j=j, drop=drop))
# Line 354  Line 416 
416  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## 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*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
418  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
419      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
420      if(missing(i))      if(missing(i))
421          x[, j] <- value          x[, j] <- value
422      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 364  Line 426 
426              x[i, ] <- value              x[i, ] <- value
427          else if(na == 3)          else if(na == 3)
428              x[i] <- value              x[i] <- value
429          else stop("Internal bug: nargs()=",na,"; please report")          else stop(gettextf("Internal bug: nargs()=%d; please report",
430                               na), domain=NA)
431      } else      } else
432          x[i,j] <- value          x[i,j] <- value
433      if(isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 383  Line 446 
446                           replDiag(x, i=i, , value=value)                           replDiag(x, i=i, , value=value)
447                   })                   })
448    
449    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",
450                                    j = "index", value = "replValue"),
451                     function(x,i,j, ..., value) replDiag(x, j=j, value=value))
452    
453    ## x[] <- value :
454    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",  setReplaceMethod("[", signature(x = "diagonalMatrix",
471                                  i = "matrix", # 2-col.matrix                                  i = "matrix", # 2-col.matrix
472                                  j = "missing", value = "replValue"),                                  j = "missing", value = "replValue"),
# Line 394  Line 478 
478                                   if(any(value != one | is.na(value))) {                                   if(any(value != one | is.na(value))) {
479                                       x@diag <- "N"                                       x@diag <- "N"
480                                       x@x <- rep.int(one, x@Dim[1])                                       x@x <- rep.int(one, x@Dim[1])
481                                   }                                   } else return(x)
482                               }                               }
483                               x@x[ii] <- value                               x@x[ii] <- value
484                               x                               x
485                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
486                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
487                               x[i] <- value                               x[i] <- value
488                               x                               x
489                           }                           }
# Line 411  Line 495 
495                       }                       }
496                   })                   })
497    
 setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  
                                 j = "index", value = "replValue"),  
                  function(x,i,j, ..., value) replDiag(x, j=j, value=value))  
498    
499    ## value = "sparseMatrix":
500  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
501                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
502                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
503                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
504  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
505                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
506                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
507                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
508  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
509                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
510                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
511                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
512    
513    ## value = "sparseVector":
514  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
515                                  value = "sparseVector"),                                  value = "sparseVector"),
516                   replDiag)                   replDiag)
# Line 442  Line 525 
525  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
526            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
527    
528  setMethod("isDiagonal", signature(object = "diagonalMatrix"),  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
529            function(object) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, upper=NA, ...) TRUE)
530  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
           function(object) TRUE)  
 setMethod("isSymmetric", signature(object = "diagonalMatrix"),  
           function(object, ...) TRUE)  
531    
532  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
533  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), function(x) .setZero(x))
534    
535  setMethod("chol", signature(x = "ddiMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),
536            function(x, pivot, ...) {            function(x, pivot, ...) {
# Line 484  Line 564 
564  ##       ---------------------  ##       ---------------------
565  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
566  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
567      n <- dimCheck(x,y)[1]      dimCheck(x,y)
568      if(x@diag != "U") {      if(x@diag != "U") {
569          if(y@diag != "U") {          if(y@diag != "U") {
570              nx <- x@x * y@x              nx <- x@x * y@x
# Line 497  Line 577 
577      return(y)      return(y)
578  }  }
579    
580    ##' Boolean Algebra/Arithmetic Product of Diagonal Matrices
581    ##'  %&%
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  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
598            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
599    
600  formals(diagdiagprod) <- alist(x=, y=x)  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"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
628            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
629  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
630            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
631  setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),  
632            diagdiagprod, valueClass = "ddiMatrix")  ##' crossprod(x) := x'x
633  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),  diagprod <- function(x, y = NULL, boolArith = NA, ...) {
634            diagdiagprod, valueClass = "ddiMatrix")      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      ## x is diagonalMatrix      ## x is diagonalMatrix
656      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
657      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
658      if(dx[2] != dy[1]) stop("non-matching dimensions")  }
659      n <- dx[1]  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
660      as(if(x@diag == "U") y else x@x * y, "Matrix")  
661  }  formals(diagmatprod) <- alist(x=, y=NULL, boolArith = NA, ...=) ## FIXME boolArith
662  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  diagmatprod2 <- function(x, y=NULL, boolArith = NA, ...) {
663            diagmatprod)      ## x is diagonalMatrix
664  ## sneaky .. :      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
665  formals(diagmatprod) <- alist(x=, y=NULL)      Matrix(if(x@diag == "U") y else x@x * y)
666  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  }
667            diagmatprod)  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod2)
668    
669  diagGeprod <- 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"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
676  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
 formals(diagGeprod) <- alist(x=, y=NULL)  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  
           diagGeprod, valueClass = "dgeMatrix")  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  
           diagGeprod)  
677    
678    diagGeprodBool <- function(x, y) {
679        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    
699    ## analogous to diagmatprod() above:
700  matdiagprod <- function(x, y) {  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")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
703      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
704  }  }
705  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
           matdiagprod)  
 formals(matdiagprod) <- alist(x=, y=NULL)  
 setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  
           matdiagprod)  
706    
707  gediagprod <- function(x, y) {  gediagprod <- 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(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
712      x      x
713  }  }
714  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
715  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
716  formals(gediagprod) <- alist(x=, y=NULL)  
717  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprodBool <- function(x, y) {
718            gediagprod)      dx <- dim(x)
719  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
720            gediagprod)      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    
# Line 585  Line 775 
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    ## + '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  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
869            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            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"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
878            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            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"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
884            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            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"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
890            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            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    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
898    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
899    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
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  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  
906  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
907  ## ==> do this:            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
908  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  
909            function(x, y) as(x, "CsparseMatrix") %*% y)  setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
910              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
911    
912    }
913    
914  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
915            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
916  ## NB: this is *not* needed for Tsparse & Rsparse  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*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
920  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
921    
   
   
922  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
923            function(a, b, ...) {            function(a, b, ...) {
924                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 628  Line 927 
927            })            })
928    
929  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
930      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
931          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
932      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
933      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 647  Line 946 
946  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
947    
948  ## Use function for several signatures, in order to evade  ## Use function for several signatures, in order to evade
 ## ambiguous dispatch for "ddi", since there's also Arith(ddense., ddense.)  
949  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
950      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
951      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"      r <- callGeneric(.diag.x(e1), .diag.x(e2)) # error if not "compatible"
# Line 655  Line 953 
953      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
954                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
955      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
956          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
957              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
958                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
959              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
960                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
961                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
962                if(!is.double(r)) r <- as.double(r)
963          }          }
964          else if(is.logical(r))          else if(is.logical(r))
965              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
966          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
967                               typeof(r)), domain=NA)
968          e1@x <- r          e1@x <- r
969          .diag.2N(e1)          .diag.2N(e1)
970      }      }
971      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
972            ## e.g., m == m
973          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
974          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
975            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
976          d <- e1@Dim          d <- e1@Dim
977          n <- d[1]          n <- d[1]
978          stopifnot(length(r) == n)          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))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
981          newcl <-          newcl <-
982              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
983                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
984              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
985    
986          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
987      }      }
# Line 701  Line 999 
999              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
1000  }  }
1001    
1002  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
1003  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
1004  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
1005  ##     before conversion to sparse, since that loses "diagonality"}  ##     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" :  ## For almost everything else, diag* shall be treated "as sparse" :
1067  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 728  Line 1085 
1085  ##       Compare --> logical  ##       Compare --> logical
1086  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
1087    
1088  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
1089  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
1090  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1091    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
1092            function(e1,e2) {            function(e1,e2) {
1093                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1094                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1095                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1096                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1097                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
1098                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
1099                        e1@diag <- "N"                        e1@diag <- "N"
1100                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1101                    } else                        } ## else: result = e1  (is "U" diag)
1102                      } else {
1103                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1104                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1105                    return(e1)                        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)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
1110            })            })
1111    
1112  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
1113    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
1114            function(e1,e2) {            function(e1,e2) {
1115                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
1116                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1117                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1118                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1119                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
1120                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1121                        e2@diag <- "N"                        e2@diag <- "N"
1122                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
1123                    } else                        } ## else: result = e2  (is "U" diag)
1124                      } else {
1125                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1126                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1127                    return(e2)                        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"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
1132            })            })
1133    
1134  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
1135  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1136    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
1137            function(e1,e2) {            function(e1,e2) {
1138                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1139                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1140                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1141                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1142                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1143                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "double"
1144                        e1@diag <- "N"                    if(e1@diag == "U") {
1145                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1146                    } else                            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)                        r <- callGeneric(e1@x, e2)
1151                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1152                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
1153                }                }
1154                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
1155                  } else
1156                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1157            })            })
1158    
1159  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
1160    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
1161            function(e1,e2) {            function(e1,e2) {
1162                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
1163                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1164                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1165                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1166                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1167                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## storage.mode(E@x) <- "double"
1168                        e2@diag <- "N"                    if(e2@diag == "U") {
1169                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1170                    } else                            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)                        r <- callGeneric(e1, e2@x)
1175                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1176                    e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e2)  
1177                }                }
1178                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
1179                  } else
1180                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1181            })            })
1182    
1183  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1184  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
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) {            function(e1,e2) {
1194                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1195                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1196                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1197                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1198                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1199                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "logical"
1200                        e1@diag <- "N"                    if(e1@diag == "U") {
1201                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1202                    } else                            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)                        r <- callGeneric(e1@x, e2)
1207                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1208                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        E@x[seq_len(n)] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
                   return(e1)  
1209                }                }
1210                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    E
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n^2  
               f0 <- callGeneric(e1, 0)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1211                    } else                    } else
1212                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "d"), e2)
                   e2 <- copyClass(e2, "ldiMatrix", c("diag", "Dim", "Dimnames"))  
                   e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1213            })            })
1214    
1215  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1216  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1217    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1218            function(e1,e2) {            function(e1,e2) {
1219                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1220                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1221                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1222                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1223                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1224                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1225                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1226                        e1@diag <- "N"                        e1@diag <- "N"
1227                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1228                    } else                        } ## else: result = e1  (is "U" diag)
1229                      } else {
1230                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1231                    e1@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1232                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1233                }                }
1234                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n^2  
               f0 <- callGeneric(e1, FALSE)  
               if(all(is0(f0))) { # remain diagonal  
                   L1 <- (le <- length(e1)) == 1L  
                   if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)  
                   if(e2@diag == "U" && any((r <- callGeneric(e1, TRUE)) != 1)) {  
                       e2@diag <- "N"  
                       if(L1) r <- rep.int(r, n)  
1235                    } else                    } else
1236                        r <- callGeneric(e1, e2@x)                    callGeneric(diag2tT.u(e1,e2, "l"), e2)
                   e2@x <- if(L1) r else r[1L + (n+1L)*(0:(n-1L))]  
                   return(e2)  
               }  
               callGeneric(e1, diag2tT.u(e2,e1, "l"))  
1237            })            })
1238    
1239    
   
1240  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1241  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1242      ## ddi*:      ## ddi*:
1243      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1244                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1245      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1246                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1247      ## ldi*:      ## ldi*:
1248      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1249                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1250      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1251                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1252  }  }
1253  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1254  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1255  for(cl in diCls) {  if(FALSE) # now also contains "geMatrix"
1256      setMethod("&", signature(e1 = cl, e2 = "ANY"),  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1257                function(e1,e2) e1 & as(e2,"Matrix"))                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1258      setMethod("&", signature(e1 = "ANY", e2 = cl),  dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1259                function(e1,e2) as(e1,"Matrix") & e2)  for(DI in diCls) {
1260      for(c2 in c("denseMatrix", "Matrix")) {      dMeth <- if(extends(DI, "dMatrix"))
1261          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1262                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      else # "lMatrix", the only other kind for now
1263          setMethod("&", signature(e1 = c2, e2 = cl),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1264                    function(e1,e2) Diagonal(x = diag(e1)) & e2)      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"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1281  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1282  for(cl in diCls) {  for(cl in diCls) {
1283  setMethod("any", cl,  setMethod("any", cl,
1284            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1285                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1286                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)                else if(x@diag == "U") TRUE else any(x@x, ..., na.rm = na.rm)
1287            })            })
1288  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1289  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      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,  setMethod("sum", cl,
1303            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 960  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) {            function(object) {
1344                d <- dim(object)                d <- dim(object)
# Line 975  Line 1354 
1354                    invisible(object)                    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.2226  
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