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 2840, Fri Oct 5 22:18:37 2012 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))      if(missing(n))
10          n <- length(x)          n <- length(x)
11      else {      else {
# Line 17  Line 17 
17          new("ddiMatrix", Dim = c(n,n), diag = "U")          new("ddiMatrix", Dim = c(n,n), diag = "U")
18      else {      else {
19          lx <- length(x)          lx <- length(x)
20          stopifnot(lx == 1 || lx == n) # but keep 'x' short for now          lx.1 <- lx == 1L
21            stopifnot(lx.1 || lx == n) # but keep 'x' short for now
22          if(is.logical(x))          if(is.logical(x))
23              cl <- "ldiMatrix"              cl <- "ldiMatrix"
24          else if(is.numeric(x)) {          else if(is.numeric(x)) {
# Line 27  Line 28 
28          else if(is.complex(x)) {          else if(is.complex(x)) {
29              cl <- "zdiMatrix"  # will not yet work              cl <- "zdiMatrix"  # will not yet work
30          } else stop("'x' has invalid data type")          } else stop("'x' has invalid data type")
31            if(lx.1 && !is.na(x) && x == 1) # cheap check for uni-diagonal..
32                new(cl, Dim = c(n,n), diag = "U")
33            else
34          new(cl, Dim = c(n,n), diag = "N",          new(cl, Dim = c(n,n), diag = "N",
35              x = if(lx == 1) rep.int(x,n) else x)                  x = if(lx.1) rep.int(x,n) else x)
36      }      }
37  }  }
38    
39  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {  .sparseDiagonal <- function(n, x = 1, uplo = "U",
40                                shape = if(missing(cols)) "t" else "g",
41                                unitri, kind,
42                                cols = if(n) 0:(n - 1L) else integer(0))
43    {
44      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
45      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!(mcols <- missing(cols)))
46      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
47      stopifnot(is.character(shape), nchar(shape) == 1,      m <- length(cols)
48                any(shape == c("t","s","g"))) # triangular / symmetric / general      if(missing(kind))
49      kind <-      kind <-
50          if(is.double(x)) "d"          if(is.double(x)) "d"
51          else if(is.logical(x)) "l"          else if(is.logical(x)) "l"
# Line 45  Line 53 
53              storage.mode(x) <- "double"              storage.mode(x) <- "double"
54              "d"              "d"
55          }          }
56      new(paste(kind, shape, "CMatrix", sep=''),      else stopifnot(any(kind == c("d","l","n")))
57          Dim = c(n,n), x = x, uplo = uplo,      stopifnot(is.character(shape), nchar(shape) == 1,
58          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)                any(shape == c("t","s","g"))) # triangular / symmetric / general
59        if((missing(unitri) || unitri) && shape == "t" &&
60           (mcols || cols == 0:(n-1L)) &&
61           ((any(kind == c("l", "n")) && allTrue(x)) ||
62            (    kind == "d"          && allTrue(x == 1)))) { ## uni-triangular
63            new(paste0(kind,"tCMatrix"), Dim = c(n,n),
64                       uplo = uplo, diag = "U", p = rep.int(0L, n+1L))
65        }
66        else if(kind == "n") {
67            if(shape == "g")
68                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
69            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
70                     i = cols, p = 0:m)
71        }
72        else { ## kind != "n" -- have x slot :
73            if((lx <- length(x)) == 1) x <- rep.int(x, m)
74            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
75            if(shape == "g")
76                new(paste0(kind, "gCMatrix"), Dim = c(n,m),
77                    x = x, i = cols, p = 0:m)
78            else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
79                     x = x, i = cols, p = 0:m)
80        }
81  }  }
82    
83  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
84  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
85      .sparseDiagonal(n, x, uplo, shape = "s")      .sparseDiagonal(n, x, uplo, shape = "s")
86    
87  ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
88  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
89      .sparseDiagonal(n, x, uplo, shape = "t")      .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
90    
91    
92  ### 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.
93  ### 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
94  ### 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
95  ### who posted his bdiag() function written in December 1995.  ## who posted his bdiag() function written in December 1995.
96  if(FALSE)##--- no longer used:  if(FALSE)##--- no longer used:
97  .bdiag <- function(lst) {  .bdiag <- function(lst) {
98      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
99      stopifnot(is.list(lst), length(lst) >= 1)      stopifnot(is.list(lst), length(lst) >= 1)
100      dims <- sapply(lst, dim, USE.NAMES=FALSE)      dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
101      ## make sure we had all matrices:      ## make sure we had all matrices:
102      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
103          stop("some arguments are not matrices")          stop("some arguments are not matrices")
# Line 96  Line 126 
126                     as, "TsparseMatrix")                     as, "TsparseMatrix")
127      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
128      ## else      ## else
129      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))      i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
130      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))      j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
131    
132      clss <- sapply(Tlst, class)      clss <- vapply(Tlst, class, "")
133      knds <- substr(clss, 2, 2)      typ <- substr(clss, 2, 2)
134      sym  <- knds == "s" # symmetric ones      knd <- substr(clss, 1, 1)
135      tri  <- knds == "t" # triangular ones      sym <- typ == "s" # symmetric ones
136      use.n <- any(is.n <- substr(clss,1,1) == "n")      tri <- typ == "t" # triangular ones
137      if(use.n && !(use.n <- all(is.n)))      use.n <- any(is.n <- knd == "n")
138        if(use.n && !(use.n <- all(is.n))) {
139          Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")          Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
140            knd [is.n] <- "l"
141        }
142        use.l <- !use.n && all(knd == "l")
143      if(all(sym)) { ## result should be *symmetric*      if(all(sym)) { ## result should be *symmetric*
144          uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"          uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
145          tLU <- table(uplos)# of length 1 or 2 ..          tLU <- table(uplos)# of length 1 or 2 ..
146          if(length(tLU) == 1) { ## all "U" or all "L"          if(length(tLU) == 1) { ## all "U" or all "L"
147              useU <- uplos[1] == "U"              useU <- uplos[1] == "U"
# Line 121  Line 155 
155          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
156              r <- new("nsTMatrix")              r <- new("nsTMatrix")
157          } else {          } else {
158              r <- new("dsTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
159              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
160          }          }
161          r@uplo <- if(useU) "U" else "L"          r@uplo <- if(useU) "U" else "L"
162      }      }
163      else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"      else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
164                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower                            all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
165         ){ ## *triangular* result         ){ ## *triangular* result
166    
167          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
168              r <- new("ntTMatrix")              r <- new("ntTMatrix")
169          } else {          } else {
170              r <- new("dtTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
171              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
172          }          }
173          r@uplo <- ULs[1L]          r@uplo <- ULs[1L]
# Line 144  Line 178 
178          if(use.n) { ## return nsparseMatrix :          if(use.n) { ## return nsparseMatrix :
179              r <- new("ngTMatrix")              r <- new("ngTMatrix")
180          } else {          } else {
181              r <- new("dgTMatrix")              r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
182              r@x <- unlist(lapply(Tlst, slot, "x"))              r@x <- unlist(lapply(Tlst, slot, "x"))
183          }          }
184      }      }
# Line 170  Line 204 
204  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
205      ## to triangular Tsparse      ## to triangular Tsparse
206      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
207      new(paste(kind, "tTMatrix", sep=''),      new(paste0(kind, "tTMatrix"),
208          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,          diag = from@diag, Dim = from@Dim, Dimnames = from@Dimnames,
209          uplo = uplo,          uplo = uplo,
210          x = from@x, # <- ok for diag = "U" and "N" (!)          x = from@x, # <- ok for diag = "U" and "N" (!)
# Line 181  Line 215 
215      ## to symmetric Tsparse      ## to symmetric Tsparse
216      n <- from@Dim[1]      n <- from@Dim[1]
217      i <- seq_len(n) - 1L      i <- seq_len(n) - 1L
218      new(paste(kind, "sTMatrix", sep=''),      new(paste0(kind, "sTMatrix"),
219          Dim = from@Dim, Dimnames = from@Dimnames,          Dim = from@Dim, Dimnames = from@Dimnames,
220          i = i, j = i, uplo = uplo,          i = i, j = i, uplo = uplo,
221          x = if(from@diag == "N") from@x else ## "U"-diag          x = if(from@diag == "N") from@x else ## "U"-diag
# Line 206  Line 240 
240          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
241  }  }
242    
243    ## FIXME: should not be needed {when ddi* is dsparse* etc}:
244    setMethod("is.finite", signature(x = "diagonalMatrix"),
245              function(x) is.finite(.diag2tT(x)))
246    setMethod("is.infinite", signature(x = "diagonalMatrix"),
247              function(x) is.infinite(.diag2tT(x)))
248    
249  ## In order to evade method dispatch ambiguity warnings,  ## In order to evade method dispatch ambiguity warnings,
250  ## 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
# Line 214  Line 253 
253  ## ddi*:  ## ddi*:
254  diag2tT <- function(from) .diag2tT(from, "U", "d")  diag2tT <- function(from) .diag2tT(from, "U", "d")
255  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", diag2tT)
256  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
257  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
258  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
259    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
260  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
261        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
262  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 225  Line 265 
265  ## ldi*:  ## ldi*:
266  diag2tT <- function(from) .diag2tT(from, "U", "l")  diag2tT <- function(from) .diag2tT(from, "U", "l")
267  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", diag2tT)
268  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
269  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
270  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
271    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
272  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
273        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
274  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 247  Line 288 
288  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
289  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
290      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
291      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
292      y      y
293  }  }
294    
# Line 264  Line 305 
305                mod.x <- mode(x@x)                mod.x <- mode(x@x)
306                r <- vector(mod.x, length = n^2)                r <- vector(mod.x, length = n^2)
307                if(n)                if(n)
308                    r[1 + 0:(n - 1) * (n + 1)] <-                    r[1 + 0:(n - 1L) * (n + 1)] <-
309                        if(x@diag == "U") as1(mod=mod.x) else x@x                        if(x@diag == "U") as1(mod=mod.x) else x@x
310                r                r
311            })            })
# Line 272  Line 313 
313  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:  setAs("diagonalMatrix", "generalMatrix", # prefer sparse:
314        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))        function(from) as(as(from, "CsparseMatrix"), "generalMatrix"))
315    
316    setAs("diagonalMatrix", "denseMatrix",
317          function(from) as(as(from, "CsparseMatrix"), "denseMatrix"))
318    
319  .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
320    
321  .diag.2N <- function(m) {  .diag.2N <- function(m) {
# Line 292  Line 336 
336            d <- dim(from)            d <- dim(from)
337            if(d[1] != (n <- d[2])) stop("non-square matrix")            if(d[1] != (n <- d[2])) stop("non-square matrix")
338            if(any(from[row(from) != col(from)] != 0))            if(any(from[row(from) != col(from)] != 0))
339                stop("has non-zero off-diagonal entries")                stop("matrix with non-zero off-diagonals cannot be coerced to diagonalMatrix")
340            x <- diag(from)            x <- diag(from)
341            if(is.logical(x)) {            if(is.logical(x)) {
342                cl <- "ldiMatrix"                cl <- "ldiMatrix"
343                uni <- all(x)                uni <- allTrue(x) ## uni := {is it unit-diagonal ?}
344            } else {            } else {
345                cl <- "ddiMatrix"                cl <- "ddiMatrix"
346                uni <- all(x == 1)                uni <- allTrue(x == 1)
347                storage.mode(x) <- "double"                storage.mode(x) <- "double"
348            } ## TODO: complex            } ## TODO: complex
349            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 360 
360            x <- diag(from)            x <- diag(from)
361            if(is.logical(x)) {            if(is.logical(x)) {
362                cl <- "ldiMatrix"                cl <- "ldiMatrix"
363                uni <- all(x)                uni <- allTrue(x)
364            } else {            } else {
365                cl <- "ddiMatrix"                cl <- "ddiMatrix"
366                uni <- all(x == 1)                uni <- allTrue(x == 1)
367                storage.mode(x) <- "double"                storage.mode(x) <- "double"
368            }            } ## TODO: complex
369            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",            new(cl, Dim = c(n,n), diag = if(uni) "U" else "N",
370                x = if(uni) x[FALSE] else x)                x = if(uni) x[FALSE] else x)
371        })        })
# Line 331  Line 375 
375            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
376    
377  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
378      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
379      x <- if(missing(i))      x <- if(missing(i))
380          x[, j, drop=drop]          x[, j, drop=drop]
381      else if(missing(j))      else if(missing(j))
382          x[i, , drop=drop]          if(nargs() == 4) x[i, , drop=drop] else x[i, drop=drop]
383      else      else
384          x[i,j, drop=drop]          x[i,j, drop=drop]
385      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x      if(isS4(x) && isDiagonal(x)) as(x, "diagonalMatrix") else x
# Line 345  Line 389 
389                           j = "index", drop = "logical"), subDiag)                           j = "index", drop = "logical"), subDiag)
390  setMethod("[", signature(x = "diagonalMatrix", i = "index",  setMethod("[", signature(x = "diagonalMatrix", i = "index",
391                          j = "missing", drop = "logical"),                          j = "missing", drop = "logical"),
392            function(x, i, j, ..., drop) subDiag(x, i=i, drop=drop))            function(x, i, j, ..., drop) {
393                  na <- nargs()
394                  Matrix.msg("diag[i,m,l] : nargs()=", na, .M.level = 2)
395                  if(na == 4)
396                       subDiag(x, i=i, , drop=drop)
397                  else subDiag(x, i=i,   drop=drop)
398              })
399  setMethod("[", signature(x = "diagonalMatrix", i = "missing",  setMethod("[", signature(x = "diagonalMatrix", i = "missing",
400                           j = "index", drop = "logical"),                           j = "index", drop = "logical"),
401            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 405 
405  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch  ## FIXME: this now fails because the "denseMatrix" methods come first in dispatch
406  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*  ## Only(?) current bug:  x[i] <- value  is wrong when  i is *vector*
407  replDiag <- function(x, i, j, ..., value) {  replDiag <- function(x, i, j, ..., value) {
408      x <- as(x, "sparseMatrix")      x <- as(x, "CsparseMatrix")# was "Tsparse.." till 2012-07
409      if(missing(i))      if(missing(i))
410          x[, j] <- value          x[, j] <- value
411      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v      else if(missing(j)) { ##  x[i , ] <- v  *OR*   x[i] <- v
# Line 395  Line 445 
445                                   if(any(value != one | is.na(value))) {                                   if(any(value != one | is.na(value))) {
446                                       x@diag <- "N"                                       x@diag <- "N"
447                                       x@x <- rep.int(one, x@Dim[1])                                       x@x <- rep.int(one, x@Dim[1])
448                                   }                                   } else return(x)
449                               }                               }
450                               x@x[ii] <- value                               x@x[ii] <- value
451                               x                               x
452                           } else { ## no longer diagonal, but remain sparse:                           } else { ## no longer diagonal, but remain sparse:
453                               x <- as(x, "sparseMatrix")                               x <- as(x, "TsparseMatrix")
454                               x[i] <- value                               x[i] <- value
455                               x                               x
456                           }                           }
# Line 417  Line 467 
467                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))
468    
469  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
470                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
471                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
472                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
473  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
474                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
475                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
476                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
477  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
478                                  value = "scarceMatrix"),                                  value = "sparseMatrix"),
479                   function (x, i, j, ..., value)                   function (x, i, j, ..., value)
480                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))                   callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
481    
# Line 485  Line 535 
535  ##       ---------------------  ##       ---------------------
536  ## Note that "ldi" logical are treated as numeric  ## Note that "ldi" logical are treated as numeric
537  diagdiagprod <- function(x, y) {  diagdiagprod <- function(x, y) {
538      n <- dimCheck(x,y)[1]      dimCheck(x,y)
539      if(x@diag != "U") {      if(x@diag != "U") {
540          if(y@diag != "U") {          if(y@diag != "U") {
541              nx <- x@x * y@x              nx <- x@x * y@x
# Line 517  Line 567 
567      dx <- dim(x)      dx <- dim(x)
568      dy <- dim(y)      dy <- dim(y)
569      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
     n <- dx[1]  
570      as(if(x@diag == "U") y else x@x * y, "Matrix")      as(if(x@diag == "U") y else x@x * y, "Matrix")
571  }  }
572  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
# Line 586  Line 635 
635  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
636  ##           })  ##           })
637    
638  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  Cspdiagprod <- function(x, y) {
639  ##        function(x, y = NULL) {      dx <- dim(x <- .Call(Csparse_diagU2N, x))
640  ##           })      dy <- dim(y)
641        if(dx[2] != dy[1]) stop("non-matching dimensions")
642        if(y@diag == "N") {
643            if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
644                x <- as(x, "generalMatrix")
645            ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
646            x@x <- x@x * y@x[ind]
647        }
648        if(is(x, "compMatrix") && length(xf <- x@factors)) {
649            ## instead of dropping all factors, be smart about some
650            ## TODO ......
651            x@factors <- list()
652        }
653        x
654    }
655    
656    diagCspprod <- function(x, y) {
657        dx <- dim(x)
658        dy <- dim(y <- .Call(Csparse_diagU2N, y))
659        if(dx[2] != dy[1]) stop("non-matching dimensions")
660        if(x@diag == "N") {
661            if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
662                y <- as(y, "generalMatrix")
663            y@x <- y@x * x@x[y@i + 1L]
664        }
665        if(is(y, "compMatrix") && length(yf <- y@factors)) {
666            ## instead of dropping all factors, be smart about some
667            ## TODO
668            keep <- character()
669            if(iLU <- names(yf) == "LU") {
670                ## TODO keep <- "LU"
671            }
672            y@factors <- yf[keep]
673        }
674        y
675    }
676    
677    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
678              function(x, y = NULL) diagCspprod(x, y))
679    
680  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
681            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))
682    
683    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
684    ##  x'y == (y'x)'
685    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
686              function(x, y = NULL) t(diagCspprod(y, x)))
687    
688  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
689            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))
690    
691    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
692              function(x, y = NULL) diagCspprod(x, t(y)))
693    
694  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
695            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))
696    
697    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
698              function(x, y = NULL) Cspdiagprod(x, y))
699    
700  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
701            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))
702    
703    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
704              function(x, y) diagCspprod(x, y))
705    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
706  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
707            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))
708    
709  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
710            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))
711  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  
 ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  
 ## ==> do this:  
 setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  
           function(x, y) as(x, "CsparseMatrix") %*% y)  
712  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
713            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y))
714  ## NB: this is *not* needed for Tsparse & Rsparse  
715  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
716  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
717    
   
   
718  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
719            function(a, b, ...) {            function(a, b, ...) {
720                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 629  Line 723 
723            })            })
724    
725  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
726      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
727          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
728      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
729      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 648  Line 742 
742  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
743    
744  ## 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.)  
745  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
746      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
747      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 749 
749      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
750                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
751      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
752          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
753              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
754                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
755              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
756                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
757                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
758                if(!is.double(r)) r <- as.double(r)
759          }          }
760          else if(is.logical(r))          else if(is.logical(r))
761              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
# Line 670  Line 764 
764          .diag.2N(e1)          .diag.2N(e1)
765      }      }
766      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
767            ## e.g., m == m
768          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
769          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
770            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
771          d <- e1@Dim          d <- e1@Dim
772          n <- d[1]          n <- d[1]
773          stopifnot(length(r) == n)          stopifnot(length(r) == n)
774            if(isNum && !is.double(r)) r <- as.double(r)
775          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
776          newcl <-          newcl <-
777              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
778                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
779              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
780    
781          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
782      }      }
# Line 702  Line 794 
794              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
795  }  }
796    
797  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
798  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
799  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
800  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
801    diagOtri <- function(e1,e2) {
802        ## result must be triangular
803        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
804        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
805        e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE
806        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
807        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
808            diag(e2) <- r
809            ## check what happens "in the triangle"
810            e2.2 <- if(.n2) 2 else TRUE
811            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
812                n <- dim(e2)[1L]
813                it <- indTri(n, upper = (e2@uplo == "U"))
814                e2[it] <- callGeneric(e1.0, e2[it])
815            }
816            e2
817        }
818        else { ## result not triangular ---> general
819            rr <- as(e2, "generalMatrix")
820            diag(rr) <- r
821            rr
822        }
823    }
824    
825    
826    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
827              diagOtri)
828    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
829    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
830    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
831              function(e1,e2)
832          { ## this must only trigger for *dense* e1
833              switch(.Generic,
834                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
835                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
836                     "*" = {
837                         n <- e2@Dim[1L]
838                         d2 <- if(e2@diag == "U") { # unit-diagonal
839                             d <- rep.int(as1(e2@x), n)
840                             e2@x <- d
841                             e2@diag <- "N"
842                             d
843                         } else e2@x
844                         e2@x <- diag(e1) * d2
845                         e2
846                     },
847                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
848                         e1 ^ as(e2, "denseMatrix")
849                     },
850                     ## otherwise:
851                     callGeneric(e1, diag2Tsmart(e2,e1)))
852    })
853    
854    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
855    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
856              .Cmp.swap)
857    ## '&' and "|'  are commutative:
858    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
859              function(e1,e2) callGeneric(e2, e1))
860    
861  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
862  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 729  Line 880 
880  ##       Compare --> logical  ##       Compare --> logical
881  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
882    
883  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
884  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
885  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
886    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
887            function(e1,e2) {            function(e1,e2) {
888                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
889                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
890                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
891                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
892                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
893                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
894                        e1@diag <- "N"                        e1@diag <- "N"
895                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
896                    } else                        } ## else: result = e1  (is "U" diag)
897                      } else {
898                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
899                    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
900                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
901                }                }
902                      e1
903                  } else
904                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
905            })            })
906    
907  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
908    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
909            function(e1,e2) {            function(e1,e2) {
910                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
911                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
912                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
913                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
914                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
915                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
916                        e2@diag <- "N"                        e2@diag <- "N"
917                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
918                    } else                        } ## else: result = e2  (is "U" diag)
919                      } else {
920                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
921                    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
922                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
923                }                }
924                      e2
925                  } else
926                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
927            })            })
928    
929  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
930  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
931    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
932            function(e1,e2) {            function(e1,e2) {
933                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
934                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
935                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
936                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
937                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
938                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ddiMatrix", check=FALSE)
939                        e1@diag <- "N"                    ## storage.mode(E@x) <- "double"
940                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
941                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
942                              E@diag <- "N"
943                              E@x[seq_len(n)] <- r # possibly recycling r
944                          } ## else: result = E  (is "U" diag)
945                      } else {
946                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
947                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
948                    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)  
949                }                }
950                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
951                  } else
952                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
953            })            })
954    
955  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
956    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
957            function(e1,e2) {            function(e1,e2) {
958                n <- e2@Dim[1]; nsq <- n^2                n <- e2@Dim[1]
959                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
960                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
961                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
962                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
963                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## E <- copyClass(e2, "ddiMatrix", check=FALSE)
964                        e2@diag <- "N"                    ## storage.mode(E@x) <- "double"
965                        if(L1) r <- rep.int(r, n)                    if(e2@diag == "U") {
966                    } else                        if(any((r <- callGeneric(e1, 1)) != 1)) {
967                              E@diag <- "N"
968                              E@x[seq_len(n)] <- r # possibly recycling r
969                          } ## else: result = E  (is "U" diag)
970                      } else {
971                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
972                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
973                    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)  
974                }                }
975                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
976                  } else
977                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
978            })            })
979    
980  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
981  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
982    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
983    if(FALSE) {
984        selectMethod("<", c("numeric","lMatrix"))# Compare
985        selectMethod("&", c("numeric","lMatrix"))# Logic
986    } ## so we don't need to define a method here :
987    
988    for(arg2 in c("numeric","logical"))
989    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
990            function(e1,e2) {            function(e1,e2) {
991                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
992                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
993                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
994                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
995                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE
996                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## E <- copyClass(e1, "ldiMatrix", check=FALSE)
997                        e1@diag <- "N"                    ## storage.mode(E@x) <- "logical"
998                        if(L1) r <- rep.int(r, n)                    if(e1@diag == "U") {
999                    } else                        if(any((r <- callGeneric(1, e2)) != 1)) {
1000                              E@diag <- "N"
1001                              E@x[seq_len(n)] <- r # possibly recycling r
1002                          } ## else: result = E  (is "U" diag)
1003                      } else {
1004                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1005                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1006                    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)  
1007                }                }
1008                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)  
1009                    } else                    } else
1010                        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"))  
1011            })            })
1012    
1013  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1014  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1015    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1016            function(e1,e2) {            function(e1,e2) {
1017                n <- e1@Dim[1]; nsq <- n^2                n <- e1@Dim[1]
1018                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1019                if(all(is0(f0))) { # remain diagonal                if(all(is0(f0))) { # remain diagonal
1020                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1021                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1022                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1023                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1024                        e1@diag <- "N"                        e1@diag <- "N"
1025                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1026                    } else                        } ## else: result = e1  (is "U" diag)
1027                      } else {
1028                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1029                    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
1030                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1031                }                }
1032                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)  
1033                    } else                    } else
1034                        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"))  
1035            })            })
1036    
1037    
   
1038  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1039  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1040      ## ddi*:      ## ddi*:
1041      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1042                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1043      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1044                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1045      ## ldi*:      ## ldi*:
1046      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1047                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1048      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1049                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1050  }  }
1051  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1052  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1053  for(cl in diCls) {  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1054      setMethod("&", signature(e1 = cl, e2 = "ANY"),                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1055                function(e1,e2) e1 & as(e2,"Matrix"))  for(DI in diCls) {
1056      setMethod("&", signature(e1 = "ANY", e2 = cl),      dMeth <- if(extends(DI, "dMatrix"))
1057                function(e1,e2) as(e1,"Matrix") & e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1058      for(c2 in c("denseMatrix", "Matrix")) {      else # "lMatrix", the only other kind for now
1059          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1060                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      for(c2 in c(dense.subCl, "Matrix")) {
1061          setMethod("&", signature(e1 = c2, e2 = cl),          for(Fun in c("*", "&")) {
1062                    function(e1,e2) Diagonal(x = diag(e1)) & e2)              setMethod(Fun, signature(e1 = DI, e2 = c2),
1063                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1064                setMethod(Fun, signature(e1 = c2, e2 = DI),
1065                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1066            }
1067            setMethod("^", signature(e1 = c2, e2 = DI),
1068                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1069            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1070                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1071      }      }
1072  }  }
1073    
1074    
1075  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1076  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1077  for(cl in diCls) {  for(cl in diCls) {
1078  setMethod("any", cl,  setMethod("any", cl,
1079            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1080                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1081                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)
1082            })            })
1083  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1084  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1085        if(n >= 2) FALSE
1086        else if(n == 0 || x@diag == "U") TRUE
1087        else all(x@x, ..., na.rm = na.rm)
1088    })
1089    setMethod("prod", cl, function (x, ..., na.rm) {
1090        n <- x@Dim[1]
1091        if(n >= 2) 0
1092        else if(n == 0 || x@diag == "U") 1
1093        else ## n == 1, diag = "N" :
1094            prod(x@x, ..., na.rm = na.rm)
1095    })
1096    
1097  setMethod("sum", cl,  setMethod("sum", cl,
1098            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 979  Line 1149 
1149                    invisible(object)                    invisible(object)
1150                }                }
1151            })            })
1152    
1153    rm(dense.subCl, diCls)# not used elsewhere
1154    
1155    setMethod("summary", signature(object = "diagonalMatrix"),
1156              function(object, ...) {
1157                  d <- dim(object)
1158                  r <- summary(object@x, ...)
1159                  attr(r, "header") <-
1160                      sprintf('%d x %d diagonal Matrix of class "%s"',
1161                              d[1], d[2], class(object))
1162                  ## use ole' S3 technology for such a simple case
1163                  class(r) <- c("diagSummary", class(r))
1164                  r
1165              })
1166    
1167    print.diagSummary <- function (x, ...) {
1168        cat(attr(x, "header"),"\n")
1169        class(x) <- class(x)[-1]
1170        print(x, ...)
1171        invisible(x)
1172    }

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

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