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 2237, Fri Jul 25 06:55:42 2008 UTC pkg/Matrix/R/diagMatrix.R revision 3072, Fri Mar 27 15:10:48 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)
# Line 244  Line 289 
289    
290  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  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 264  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 272  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 292  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 316  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 331  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 345  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 355  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 365  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 384  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 395  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 412  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 443  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 485  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 498  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)
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),  
           diagdiagprod, valueClass = "ddiMatrix")  
 setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),  
           diagdiagprod, valueClass = "ddiMatrix")  
631    
632    ##' crossprod(x) := x'x
633    diagprod <- function(x, y = NULL, boolArith = NA, ...) {
634        bool <- isTRUE(boolArith)
635        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
636        if(x@diag != "U") {
637            if(bool) {
638                nx <- x@x & y@x
639                x@x <- as.logical(nx)
640            } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
641                nx <- x@x * x@x
642                if(is.numeric(nx) && !is.numeric(x@x))
643                    x <- as(x, "dMatrix")
644                x@x <- as.numeric(nx)
645            }
646        }
647        x
648    }
649    setMethod( "crossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
650    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
651    
652    
653    ## analogous to matdiagprod() below:
654  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
655      ## 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, ...=)
662  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
           diagmatprod)  
 ## sneaky .. :  
 formals(diagmatprod) <- alist(x=, y=NULL)  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  
           diagmatprod)  
663    
664  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
665      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")  
666      if(x@diag != "U")      if(x@diag != "U")
667          y@x <- x@x * y@x          y@x <- x@x * y@x
668      y      y
669  }  }
670  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
671  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)  
672    
673    diagGeprodBool <- function(x, y) {
674        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
675        if(!is.logical(y@x)) y <- as(y, "lMatrix")
676        if(x@diag != "U")
677            y@x <- x@x & y@x
678        y
679    }
680    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
681    
682    diagGeprod2 <- function(x, y=NULL, boolArith = NA, ...) {
683        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
684        bool <- isTRUE(boolArith)
685        if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
686        if(x@diag != "U")
687            y@x <- if(bool) x@x & y@x else x@x * y@x
688        y
689    }
690    setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), diagGeprod2)
691    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), diagGeprod2)
692    
693    
694    ## analogous to diagmatprod() above:
695  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
696      dx <- dim(x)      dx <- dim(x)
697      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
698      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]))
699  }  }
700  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)  
701    
702  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
703      dx <- dim(x)      dx <- dim(x)
704      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
705      if(y@diag == "N")      if(y@diag == "N")
706          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
707      x      x
708  }  }
709  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
710  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
711  formals(gediagprod) <- alist(x=, y=NULL)  
712  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprodBool <- function(x, y) {
713            gediagprod)      dx <- dim(x)
714  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
715            gediagprod)      if(!is.logical(x@x)) x <- as(x, "lMatrix")
716        if(y@diag == "N")
717            x@x <- x@x & rep(y@x, each = dx[1])
718        x
719    }
720    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
721    
722    setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"),
723              function(x, y=NULL, boolArith = NA, ...) {
724                  dx <- dim(x)
725                  if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
726                  bool <- isTRUE(boolArith)
727                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
728                  Matrix(if(y@diag == "U") x else
729                         if(bool) x & rep(y@x, each = dx[1])
730                         else     x * rep(y@x, each = dx[1]))
731              })
732    
733    setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
734              function(x, y=NULL, boolArith = NA, ...) {
735                  dx <- dim(x)
736                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
737                  bool <- isTRUE(boolArith)
738                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
739                  Matrix(if(bool) t(rep.int(y@x, dx[2]) & x)
740                             else t(rep.int(y@x, dx[2]) * x))
741              })
742    
743    
744    gediagprod2 <- function(x, y=NULL, boolArith = NA, ...) {
745        dx <- dim(x)
746        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
747        bool <- isTRUE(boolArith)
748        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
749        if(y@diag == "N")
750            x@x <- if(bool) x@x & rep(y@x, each = dx[1])
751                   else     x@x * rep(y@x, each = dx[1])
752        x
753    }
754    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), gediagprod2)
755    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), gediagprod2)
756    
757    
758  ## crossprod {more of these}  ## crossprod {more of these}
759    
# Line 586  Line 770 
770  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
771  ##           })  ##           })
772    
773  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  ##' @param x CsparseMatrix
774  ##        function(x, y = NULL) {  ##' @param y diagonalMatrix
775  ##           })  ##' @return x %*% y
776    Cspdiagprod <- function(x, y, boolArith = NA, ...) {
777        if((m <- ncol(x)) != y@Dim[1]) stop("non-matching dimensions")
778        if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
779            x <- .Call(Csparse_diagU2N, x)
780            cx <- getClass(class(x))
781            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
782                x <- as(x, "generalMatrix")
783            ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
784            if(isTRUE(boolArith)) {
785                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
786                x@x <- x@x & y@x[x@i + 1L]
787                if(!extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
788            } else {
789                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
790                x@x <- x@x * y@x[ind]
791            }
792            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
793                ## instead of dropping all factors, be smart about some
794                ## TODO ......
795                x@factors <- list()
796            }
797        }
798        x
799    }
800    
801    ##' @param x diagonalMatrix
802    ##' @param y CsparseMatrix
803    ##' @return x %*% y
804    diagCspprod <- function(x, y, boolArith = NA, ...) {
805        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
806        if(x@diag == "N") {
807            y <- .Call(Csparse_diagU2N, y)
808            cy <- getClass(class(y))
809            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
810                y <- as(y, "generalMatrix")
811            if(isTRUE(boolArith)) {
812                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
813                y@x <- y@x & x@x[y@i + 1L]
814                if(!extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
815            } else {
816                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
817                y@x <- y@x * x@x[y@i + 1L]
818            }
819            if(.hasSlot(y, "factors") && length(y@factors)) {
820         ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
821                ## instead of dropping all factors, be smart about some
822                ## keep <- character()
823                ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
824                ##     keep <- "LU"
825                ## }
826                ## y@factors <- yf[keep]
827                y@factors <- list()
828            }
829        }
830        y
831    }
832    
833    ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
834    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
835              function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, y, boolArith=boolArith))
836    
837  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
838            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y=NULL, boolArith=NA, ...)
839                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
840    
841    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
842    ##  x'y == (y'x)'
843    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
844              function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, x, boolArith=boolArith)))
845    
846  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
847            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
848    
849    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
850              function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(y), boolArith=boolArith))
851    
852  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
853            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
854    
855    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
856              function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(x, y, boolArith=boolArith))
857    
858  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
859            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
860    
861    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
862              function(x, y) diagCspprod(x, y, boolArith=NA))
863    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
864              function(x, y) diagCspprod(x, y, boolArith=TRUE))
865    
866    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
867    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
868    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
869  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
870            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
871    
872  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
873            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
874  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  
875  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
876  ## ==> do this:            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
877  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  
878            function(x, y) as(x, "CsparseMatrix") %*% y)  setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
879              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
880    
881    }
882    
883  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
884            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
885  ## NB: this is *not* needed for Tsparse & Rsparse  setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
886              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
887    
888  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
889  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
890    
   
   
891  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
892            function(a, b, ...) {            function(a, b, ...) {
893                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 629  Line 896 
896            })            })
897    
898  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
899      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
900          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
901      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
902      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 648  Line 915 
915  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
916    
917  ## 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.)  
918  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
919      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
920      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 656  Line 922 
922      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
923                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
924      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
925          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
926              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
927                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
928              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
929                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
930                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
931                if(!is.double(r)) r <- as.double(r)
932          }          }
933          else if(is.logical(r))          else if(is.logical(r))
934              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
935          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
936                               typeof(r)), domain=NA)
937          e1@x <- r          e1@x <- r
938          .diag.2N(e1)          .diag.2N(e1)
939      }      }
940      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
941            ## e.g., m == m
942          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
943          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
944            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
945          d <- e1@Dim          d <- e1@Dim
946          n <- d[1]          n <- d[1]
947          stopifnot(length(r) == n)          stopifnot(length(r) == n)
948            if(isNum && !is.double(r)) r <- as.double(r)
949          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
950          newcl <-          newcl <-
951              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
952                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
953              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
954    
955          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
956      }      }
# Line 702  Line 968 
968              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
969  }  }
970    
971  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
972  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
973  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
974  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
975    diagOtri <- function(e1,e2) {
976        ## result must be triangular
977        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
978        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
979        e1.0 <- if(is.numeric(d1)) 0 else FALSE
980        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
981        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
982            diag(e2) <- r
983            ## check what happens "in the triangle"
984            e2.2 <- if(.n2) 2 else TRUE
985            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
986                n <- dim(e2)[1L]
987                it <- indTri(n, upper = (e2@uplo == "U"))
988                e2[it] <- callGeneric(e1.0, e2[it])
989            }
990            e2
991        }
992        else { ## result not triangular ---> general
993            rr <- as(e2, "generalMatrix")
994            diag(rr) <- r
995            rr
996        }
997    }
998    
999    
1000    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
1001              diagOtri)
1002    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
1003    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
1004    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1005              function(e1,e2)
1006          { ## this must only trigger for *dense* e1
1007              switch(.Generic,
1008                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
1009                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
1010                     "*" = {
1011                         n <- e2@Dim[1L]
1012                         d2 <- if(e2@diag == "U") { # unit-diagonal
1013                             d <- rep.int(as1(e2@x), n)
1014                             e2@x <- d
1015                             e2@diag <- "N"
1016                             d
1017                         } else e2@x
1018                         e2@x <- diag(e1) * d2
1019                         e2
1020                     },
1021                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
1022                         e1 ^ as(e2, "denseMatrix")
1023                     },
1024                     ## otherwise:
1025                     callGeneric(e1, diag2Tsmart(e2,e1)))
1026    })
1027    
1028    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
1029    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1030              .Cmp.swap)
1031    ## '&' and "|'  are commutative:
1032    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
1033              function(e1,e2) callGeneric(e2, e1))
1034    
1035  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
1036  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 729  Line 1054 
1054  ##       Compare --> logical  ##       Compare --> logical
1055  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
1056    
1057  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
1058  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
1059  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1060    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
1061            function(e1,e2) {            function(e1,e2) {
1062                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1063                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1064                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1065                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1066                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
1067                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
1068                        e1@diag <- "N"                        e1@diag <- "N"
1069                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1070                    } else                        } ## else: result = e1  (is "U" diag)
1071                      } else {
1072                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1073                    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
1074                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1075                }                }
1076                      e1
1077                  } else
1078                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
1079            })            })
1080    
1081  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
1082    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
1083            function(e1,e2) {            function(e1,e2) {
1084                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
1085                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1086                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1087                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1088                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
1089                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1090                        e2@diag <- "N"                        e2@diag <- "N"
1091                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
1092                    } else                        } ## else: result = e2  (is "U" diag)
1093                      } else {
1094                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1095                    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
1096                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1097                }                }
1098                      e2
1099                  } else
1100                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
1101            })            })
1102    
1103  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
1104  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1105    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
1106            function(e1,e2) {            function(e1,e2) {
1107                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1108                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1109                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1110                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1111                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1112                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "double"
1113                        e1@diag <- "N"                    if(e1@diag == "U") {
1114                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1115                    } else                            E@diag <- "N"
1116                              E@x[seq_len(n)] <- r # possibly recycling r
1117                          } ## else: result = E  (is "U" diag)
1118                      } else {
1119                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1120                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1121                    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)  
1122                }                }
1123                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
1124                  } else
1125                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1126            })            })
1127    
1128  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
1129    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
1130            function(e1,e2) {            function(e1,e2) {
1131                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
1132                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1133                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1134                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1135                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1136                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## storage.mode(E@x) <- "double"
1137                        e2@diag <- "N"                    if(e2@diag == "U") {
1138                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1139                    } else                            E@diag <- "N"
1140                              E@x[seq_len(n)] <- r # possibly recycling r
1141                          } ## else: result = E  (is "U" diag)
1142                      } else {
1143                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1144                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1145                    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)  
1146                }                }
1147                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
1148                  } else
1149                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1150            })            })
1151    
1152  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1153  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
1154    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1155    if(FALSE) {
1156        selectMethod("<", c("numeric","lMatrix"))# Compare
1157        selectMethod("&", c("numeric","lMatrix"))# Logic
1158    } ## so we don't need to define a method here :
1159    
1160    for(arg2 in c("numeric","logical"))
1161    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1162            function(e1,e2) {            function(e1,e2) {
1163                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1164                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1165                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1166                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1167                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1168                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "logical"
1169                        e1@diag <- "N"                    if(e1@diag == "U") {
1170                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1171                    } else                            E@diag <- "N"
1172                              E@x[seq_len(n)] <- r # possibly recycling r
1173                          } ## else: result = E  (is "U" diag)
1174                      } else {
1175                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1176                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1177                    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)  
1178                }                }
1179                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)  
1180                    } else                    } else
1181                        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"))  
1182            })            })
1183    
1184  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1185  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1186    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1187            function(e1,e2) {            function(e1,e2) {
1188                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1189                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1190                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1191                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1192                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1193                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1194                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1195                        e1@diag <- "N"                        e1@diag <- "N"
1196                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1197                    } else                        } ## else: result = e1  (is "U" diag)
1198                      } else {
1199                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1200                    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
1201                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1202                }                }
1203                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)  
1204                    } else                    } else
1205                        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"))  
1206            })            })
1207    
1208    
   
1209  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1210  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1211      ## ddi*:      ## ddi*:
1212      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1213                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1214      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1215                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1216      ## ldi*:      ## ldi*:
1217      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1218                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1219      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1220                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1221  }  }
1222  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1223  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1224  for(cl in diCls) {  if(FALSE) # now also contains "geMatrix"
1225      setMethod("&", signature(e1 = cl, e2 = "ANY"),  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1226                function(e1,e2) e1 & as(e2,"Matrix"))                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1227      setMethod("&", signature(e1 = "ANY", e2 = cl),  dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1228                function(e1,e2) as(e1,"Matrix") & e2)  for(DI in diCls) {
1229      for(c2 in c("denseMatrix", "Matrix")) {      dMeth <- if(extends(DI, "dMatrix"))
1230          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1231                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      else # "lMatrix", the only other kind for now
1232          setMethod("&", signature(e1 = c2, e2 = cl),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1233                    function(e1,e2) Diagonal(x = diag(e1)) & e2)      for(c2 in c(dense.subCl, "Matrix")) {
1234            for(Fun in c("*", "&")) {
1235                setMethod(Fun, signature(e1 = DI, e2 = c2),
1236                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1237                setMethod(Fun, signature(e1 = c2, e2 = DI),
1238                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1239            }
1240            setMethod("^", signature(e1 = c2, e2 = DI),
1241                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1242            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1243                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1244      }      }
1245  }  }
1246    
1247    ## Group methods "Math", "Math2" in                     --> ./Math.R
1248    
1249  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1250  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1251  for(cl in diCls) {  for(cl in diCls) {
1252  setMethod("any", cl,  setMethod("any", cl,
1253            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1254                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1255                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)
1256            })            })
1257  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1258  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1259        if(n >= 2) FALSE
1260        else if(n == 0 || x@diag == "U") TRUE
1261        else all(x@x, ..., na.rm = na.rm)
1262    })
1263    setMethod("prod", cl, function (x, ..., na.rm) {
1264        n <- x@Dim[1]
1265        if(n >= 2) 0
1266        else if(n == 0 || x@diag == "U") 1
1267        else ## n == 1, diag = "N" :
1268            prod(x@x, ..., na.rm = na.rm)
1269    })
1270    
1271  setMethod("sum", cl,  setMethod("sum", cl,
1272            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 979  Line 1323 
1323                    invisible(object)                    invisible(object)
1324                }                }
1325            })            })
1326    
1327    rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1328       dense.subCl, diCls)# not used elsewhere
1329    
1330    setMethod("summary", signature(object = "diagonalMatrix"),
1331              function(object, ...) {
1332                  d <- dim(object)
1333                  r <- summary(object@x, ...)
1334                  attr(r, "header") <-
1335                      sprintf('%d x %d diagonal Matrix of class "%s"',
1336                              d[1], d[2], class(object))
1337                  ## use ole' S3 technology for such a simple case
1338                  class(r) <- c("diagSummary", class(r))
1339                  r
1340              })
1341    
1342    print.diagSummary <- function (x, ...) {
1343        cat(attr(x, "header"),"\n")
1344        class(x) <- class(x)[-1]
1345        print(x, ...)
1346        invisible(x)
1347    }

Legend:
Removed from v.2237  
changed lines
  Added in v.3072

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