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 2197, Mon Jun 2 14:34:42 2008 UTC pkg/Matrix/R/diagMatrix.R revision 3069, Thu Mar 26 10:00:49 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  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = 1, uplo = "U",
38  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {                              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      cls <-      m <- length(cols)
46          if(is.double(x)) "dsCMatrix"      if(missing(kind))
47          else if(is.logical(x)) "lsCMatrix"          kind <-
48                if(is.double(x)) "d"
49                else if(is.logical(x)) "l"
50          else { ## for now          else { ## for now
51              storage.mode(x) <- "double"              storage.mode(x) <- "double"
52              "dsCMatrix"                  "d"
53                }
54        else stopifnot(any(kind == c("d","l","n")))
55        stopifnot(is.character(shape), nchar(shape) == 1,
56                  any(shape == c("t","s","g"))) # triangular / symmetric / general
57        if((missing(unitri) || unitri) && shape == "t" &&
58           (mcols || cols == 0:(n-1L)) &&
59           ((any(kind == c("l", "n")) && allTrue(x)) ||
60            (    kind == "d"          && allTrue(x == 1)))) { ## uni-triangular
61            new(paste0(kind,"tCMatrix"), Dim = c(n,n),
62                       uplo = uplo, diag = "U", p = rep.int(0L, n+1L))
63        }
64        else if(kind == "n") {
65            if(shape == "g")
66                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
67            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
68                     i = cols, p = 0:m)
69        }
70        else { ## kind != "n" -- have x slot :
71            if((lx <- length(x)) == 1) x <- rep.int(x, m)
72            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
73            if(shape == "g")
74                new(paste0(kind, "gCMatrix"), Dim = c(n,m),
75                    x = x, i = cols, p = 0:m)
76            else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
77                     x = x, i = cols, p = 0:m)
78          }          }
     new(cls, Dim = c(n,n), x = x, uplo = uplo,  
         i = if(n) 0:(n - 1L) else integer(0), p = 0:n)  
79  }  }
80    
81  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
82  ### Bert's code built on a post by Andy Liaw who most probably was influenced  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
83  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002      .sparseDiagonal(n, x, uplo, shape = "s")
 ### who posted his bdiag() function written in December 1995.  
84    
85  bdiag <- function(...) {  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
86      if(nargs() == 0) return(new("dgCMatrix"))  .trDiagonal <- function(n, x = 1, uplo = "U", unitri=TRUE)
87      ## else :      .sparseDiagonal(n, x, uplo, shape = "t", unitri=unitri)
88      mlist <- if (nargs() == 1) as.list(...) else list(...)  
89      dims <- sapply(mlist, dim)  
90    ## This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
91    ## Bert's code built on a post by Andy Liaw who most probably was influenced
92    ## by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
93    ## who posted his bdiag() function written in December 1995.
94    if(FALSE)##--- no longer used:
95    .bdiag <- function(lst) {
96        ## block-diagonal matrix [a dgTMatrix] from list of matrices
97        stopifnot(is.list(lst), length(lst) >= 1)
98        dims <- vapply(lst, dim, 1L, USE.NAMES=FALSE)
99      ## make sure we had all matrices:      ## 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")
102      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
103                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
104      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
105        r@Dim <- as.integer(csdim[nrow(csdim),])
106      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
107      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
108          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
109          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
110              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
111          else ## square matrix          else ## square matrix
112              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
113        }
114        r
115    }
116    ## expand(<mer>) needed something like bdiag() for lower-triangular
117    ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
118    ##  implementation for those; now extended and generalized:
119    .bdiag <- function(lst) {
120        ## block-diagonal matrix [a dgTMatrix] from list of matrices
121        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123        Tlst <- lapply(lapply(lst, as_Csp2), # includes "diagU2N"
124                       as, "TsparseMatrix")
125        if(nl == 1) return(Tlst[[1]])
126        ## else
127        i_off <- c(0L, cumsum(vapply(Tlst, nrow, 1L)))
128        j_off <- c(0L, cumsum(vapply(Tlst, ncol, 1L)))
129    
130        clss <- vapply(Tlst, class, "")
131        typ <- substr(clss, 2, 2)
132        knd <- substr(clss, 1, 1)
133        sym <- typ == "s" # symmetric ones
134        tri <- typ == "t" # triangular ones
135        use.n <- any(is.n <- knd == "n")
136        if(use.n && !(use.n <- all(is.n))) {
137            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
138            knd [is.n] <- "l"
139        }
140        use.l <- !use.n && all(knd == "l")
141        if(all(sym)) { ## result should be *symmetric*
142            uplos <- vapply(Tlst, slot, ".", "uplo") ## either "U" or "L"
143            tLU <- table(uplos)# of length 1 or 2 ..
144            if(length(tLU) == 1) { ## all "U" or all "L"
145                useU <- uplos[1] == "U"
146            } else { ## length(tLU) == 2, counting "L" and "U"
147                useU <- diff(tLU) >= 0
148                if(useU && (hasL <- tLU[1] > 0))
149                    Tlst[hasL] <- lapply(Tlst[hasL], t)
150                else if(!useU && (hasU <- tLU[2] > 0))
151                    Tlst[hasU] <- lapply(Tlst[hasU], t)
152      }      }
153      ## slightly debatable if we really should return Csparse.. :          if(use.n) { ## return nsparseMatrix :
154      as(ret, "CsparseMatrix")              r <- new("nsTMatrix")
155            } else {
156                r <- new(paste0(if(use.l) "l" else "d", "sTMatrix"))
157                r@x <- unlist(lapply(Tlst, slot, "x"))
158            }
159            r@uplo <- if(useU) "U" else "L"
160  }  }
161        else if(all(tri) && { ULs <- vapply(Tlst, slot, ".", "uplo")##  "U" or "L"
162                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
163           ){ ## *triangular* result
164    
165            if(use.n) { ## return nsparseMatrix :
166                r <- new("ntTMatrix")
167            } else {
168                r <- new(paste0(if(use.l) "l" else "d", "tTMatrix"))
169                r@x <- unlist(lapply(Tlst, slot, "x"))
170            }
171            r@uplo <- ULs[1L]
172        }
173        else {
174            if(any(sym))
175                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
176            if(use.n) { ## return nsparseMatrix :
177                r <- new("ngTMatrix")
178            } else {
179                r <- new(paste0(if(use.l) "l" else "d", "gTMatrix"))
180                r@x <- unlist(lapply(Tlst, slot, "x"))
181            }
182        }
183        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
184        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
185        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
186        r
187    }
188    
189    bdiag <- function(...) {
190        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
191        if(nA == 1 && !is.list(...))
192            return(as(..., "CsparseMatrix"))
193        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
194        if(length(alis) == 1)
195            return(as(alis[[1]], "CsparseMatrix"))
196    
197        ## else : two or more arguments
198        as(.bdiag(alis), "CsparseMatrix")
199    }
200    
201    setMethod("tril", "diagonalMatrix", function(x, k = 0, ...)
202        if(k >= 0) x else .setZero(x, paste0(.M.kind(x), "tCMatrix")))
203    
204    setMethod("triu", "diagonalMatrix", function(x, k = 0, ...)
205        if(k <= 0) x else  .setZero(x, paste0(.M.kind(x), "tCMatrix")))
206    
207    
208    
209  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .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 91  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 100  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 116  Line 247 
247          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)          .diag2tT(d, uplo = if(extends(clx,"triangularMatrix")) x@uplo else "U", kind)
248  }  }
249    
250    ## FIXME: should not be needed {when ddi* is dsparse* etc}:
251    setMethod("is.finite", signature(x = "diagonalMatrix"),
252              function(x) is.finite(.diag2tT(x)))
253    setMethod("is.infinite", signature(x = "diagonalMatrix"),
254              function(x) is.infinite(.diag2tT(x)))
255    
256  ## In order to evade method dispatch ambiguity warnings,  ## In order to evade method dispatch ambiguity warnings,
257  ## and because we can save a .M.kind() call, we use this explicit  ## and because we can save a .M.kind() call, we use this explicit
258  ## "hack"  instead of signature  x = "diagonalMatrix" :  ## "hack"  instead of signature  x = "diagonalMatrix" :
259  ##  ##
260  ## ddi*:  ## ddi*:
261  diag2tT <- function(from) .diag2tT(from, "U", "d")  di2tT <- function(from) .diag2tT(from, "U", "d")
262  setAs("ddiMatrix", "triangularMatrix", diag2tT)  setAs("ddiMatrix", "triangularMatrix", di2tT)
263  setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", di2tT)
264  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
265  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", di2tT)
266    setAs("ddiMatrix", "dsparseMatrix", di2tT)
267  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
268        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
269  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "d"))
       function(from) .diag2sT(from, "U", "d"))  
270  ##  ##
271  ## ldi*:  ## ldi*:
272  diag2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
273  setAs("ldiMatrix", "triangularMatrix", diag2tT)  setAs("ldiMatrix", "triangularMatrix", di2tT)
274  setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", di2tT)
275  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
276  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", di2tT)
277    setAs("ldiMatrix", "lsparseMatrix", di2tT)
278  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
279        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
282    
283  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
           n <- from@Dim[1]  
285            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            new("ntTMatrix", i = i, j = i, diag = from@diag,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
288        })        })
289    
290    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
291    
292  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
294      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
295      if (n > 0) y[1L + 0:(n - 1L) * (n + 1L)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
296      y      y
297  }  }
298    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
299    ## --> using diag() where possible, ==> .ddi2mat()
300    
301  setAs("diagonalMatrix", "matrix",  .diag2mat <- function(from)
302        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
303            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
304            mkDiag(if(from@diag == "U") as1(from@x) else from@x,  
305                   n = from@Dim[1])  .ddi2mat <- function(from)
306        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
307    
308    setAs("ddiMatrix", "matrix", .ddi2mat)
309    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
310    setAs("diagonalMatrix", "matrix", .diag2mat)
311    
312  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
313            function(x, mode) {            function(x, mode) {
# Line 173  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 181  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 201  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 225  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 240  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 254  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 264  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 274  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 293  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 304  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 321  Line 495 
495                       }                       }
496                   })                   })
497    
498  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing",  
499                                  j = "index", value = "replValue"),  ## value = "sparseMatrix":
500                   function(x,i,j, ..., value) replDiag(x, j=j, value=value))  setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
501                                    value = "sparseMatrix"),
502                     function (x, i, j, ..., value)
503                     callGeneric(x=x, , j=j, value = as(value, "sparseVector")))
504    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
505                                    value = "sparseMatrix"),
506                     function (x, i, j, ..., value)
507                     callGeneric(x=x, i=i, , value = as(value, "sparseVector")))
508    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
509                                    value = "sparseMatrix"),
510                     function (x, i, j, ..., value)
511                     callGeneric(x=x, i=i, j=j, value = as(value, "sparseVector")))
512    
513    ## value = "sparseVector":
514    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "missing", j = "index",
515                                    value = "sparseVector"),
516                     replDiag)
517    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "missing",
518                                    value = "sparseVector"),
519                     replDiag)
520    setReplaceMethod("[", signature(x = "diagonalMatrix", i = "index", j = "index",
521                                    value = "sparseVector"),
522                     replDiag)
523    
524    
525  setMethod("t", signature(x = "diagonalMatrix"),  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 371  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 384  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    formals(diagdiagprod) <- alist(x=, y=NULL, ...=)
604    ##                                  -----  ---  matching the [t]crossprod generic
605  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
606            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
607  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
608            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
609    
610    ##' crossprod(x) := x'x
611    diagprod <- function(x, y=NULL, ...) {
612        if(x@diag != "U") {
613            nx <- x@x * x@x
614            if(is.numeric(nx) && !is.numeric(x@x))
615                x <- as(x, "dMatrix")
616            x@x <- as.numeric(nx)
617        }
618        x
619    }
620  setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),
621            diagdiagprod, valueClass = "ddiMatrix")            diagprod, valueClass = "ddiMatrix")
622  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),
623            diagdiagprod, valueClass = "ddiMatrix")            diagprod, valueClass = "ddiMatrix")
624    
625    
626    ## analogous to matdiagprod() below:
627  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
628      ## x is diagonalMatrix      ## x is diagonalMatrix
629      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
630      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
     n <- dx[1]  
     as(if(x@diag == "U") y else x@x * y, "Matrix")  
631  }  }
632  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
633            diagmatprod)            diagmatprod)
634  ## sneaky .. :  
635  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL, ...=)
636  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
           diagmatprod)  
637    
638  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
639      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")  
640      if(x@diag != "U")      if(x@diag != "U")
641          y@x <- x@x * y@x          y@x <- x@x * y@x
642      y      y
643  }  }
644  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
645  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
646  formals(diagGeprod) <- alist(x=, y=NULL)  formals(diagGeprod) <- alist(x=, y=NULL, ...=)
647  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),
648            diagGeprod, valueClass = "dgeMatrix")            diagGeprod, valueClass = "dgeMatrix")
649  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
650            diagGeprod)            diagGeprod)
651    
652    diagGeprodBool <- function(x, y) {
653        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
654        if(!is.logical(y@x)) y <- as(y, "lMatrix")
655        if(x@diag != "U")
656            y@x <- x@x & y@x
657        y
658    }
659    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
660    
661    
662    ## analogous to diagmatprod() above:
663  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
664      dx <- dim(x)      dx <- dim(x)
665      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
666      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]))
667  }  }
668  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
669            matdiagprod)  formals(matdiagprod) <- alist(x=, y=NULL, ...=)
670  formals(matdiagprod) <- alist(x=, y=NULL)  setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
671  setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
672            matdiagprod)            function(x, y=NULL, ...) {
673                  dx <- dim(x)
674                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
675                  Matrix(t(rep.int(y@x, dx[2]) * x))
676              })
677    
678  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
679      dx <- dim(x)      dx <- dim(x)
680      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
681      if(y@diag == "N")      if(y@diag == "N")
682          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
683      x      x
684  }  }
685  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
686  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
687  formals(gediagprod) <- alist(x=, y=NULL)  formals(gediagprod) <- alist(x=, y=NULL, ...=)
688  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),
689            gediagprod)            gediagprod)
690  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),
691            gediagprod)            gediagprod)
692    
693    gediagprodBool <- function(x, y) {
694        dx <- dim(x)
695        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
696        if(!is.logical(x@x)) x <- as(x, "lMatrix")
697        if(y@diag == "N")
698            x@x <- x@x & rep(y@x, each = dx[1])
699        x
700    }
701    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
702    
703  ## crossprod {more of these}  ## crossprod {more of these}
704    
705  ## tcrossprod --- all are not yet there: do the dense ones here:  ## tcrossprod --- all are not yet there: do the dense ones here:
# Line 472  Line 715 
715  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
716  ##           })  ##           })
717    
718  ## setMethod("tcrossprod", signature(x = "denseMatrix", y = "diagonalMatrix"),  ##' @param x CsparseMatrix
719  ##        function(x, y = NULL) {  ##' @param y diagonalMatrix
720  ##           })  ##' @return x %*% y
721    Cspdiagprod <- function(x, y, boolArith = NA) {
722        if((m <- ncol(x)) != y@Dim[1]) stop("non-matching dimensions")
723        if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
724            x <- .Call(Csparse_diagU2N, x)
725            cx <- getClass(class(x))
726            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
727                x <- as(x, "generalMatrix")
728            ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
729            if(isTRUE(boolArith)) {
730                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
731                x@x <- x@x & y@x[x@i + 1L]
732                if(!extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
733            } else {
734                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
735                x@x <- x@x * y@x[ind]
736            }
737            if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
738                ## instead of dropping all factors, be smart about some
739                ## TODO ......
740                x@factors <- list()
741            }
742        }
743        x
744    }
745    
746    ##' @param x diagonalMatrix
747    ##' @param y CsparseMatrix
748    ##' @return x %*% y
749    diagCspprod <- function(x, y, boolArith = FALSE) {
750        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
751        if(x@diag == "N") {
752            y <- .Call(Csparse_diagU2N, y)
753            cy <- getClass(class(y))
754            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
755                y <- as(y, "generalMatrix")
756            if(isTRUE(boolArith)) {
757                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
758                y@x <- y@x & x@x[y@i + 1L]
759                if(!extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
760            } else {
761                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
762                y@x <- y@x * x@x[y@i + 1L]
763            }
764            if(.hasSlot(y, "factors") && length(y@factors)) {
765         ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
766                ## instead of dropping all factors, be smart about some
767                ## keep <- character()
768                ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
769                ##     keep <- "LU"
770                ## }
771                ## y@factors <- yf[keep]
772                y@factors <- list()
773            }
774        }
775        y
776    }
777    
778    ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
779    setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
780              function(x, y, boolArith=NA) diagCspprod(x, y, boolArith=boolArith))
781    
782  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
783            function(x, y = NULL) crossprod(as(x, "sparseMatrix"), y))            function(x, y, boolArith=NA)
784                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
785    
786    ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
787    ##  x'y == (y'x)'
788    setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
789              function(x, y, boolArith=NA) t(diagCspprod(y, x, boolArith=boolArith)))
790    
791  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
792            function(x, y = NULL) crossprod(x, as(y, "sparseMatrix")))            function(x, y, boolArith=NA) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
793    
794    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
795              function(x, y, boolArith=NA) diagCspprod(x, t(y), boolArith=boolArith))
796    
797  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
798            function(x, y = NULL) tcrossprod(as(x, "sparseMatrix"), y))            function(x, y, boolArith=NA) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
799    
800    setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
801              function(x, y, boolArith=NA) Cspdiagprod(x, y, boolArith=boolArith))
802    
803  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
804            function(x, y = NULL) tcrossprod(x, as(y, "sparseMatrix")))            function(x, y, boolArith=NA) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
805    
806    setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
807              function(x, y) diagCspprod(x, y, boolArith=NA))
808    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
809              function(x, y) diagCspprod(x, y, boolArith=TRUE))
810    
811    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
812    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
813    
 ## FIXME?: In theory, this can be done *FASTER*, in some cases, via tapply1()  
814  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
815            function(x, y) as(x, "sparseMatrix") %*% y)            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
816    
817  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
818            function(x, y) x %*% as(y, "sparseMatrix"))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
819  ## NB: The previous is *not* triggering for  "ddi" o "dgC" (= distance 3)  
820  ##     since there's a "ddense" o "Csparse" at dist. 2 => triggers first.  setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
821  ## ==> do this:            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
822  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  
823            function(x, y) as(x, "CsparseMatrix") %*% y)  setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
824              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
825    
826    }
827    
828  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
829            function(x, y) x %*% as(y, "CsparseMatrix"))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
830  ## NB: this is *not* needed for Tsparse & Rsparse  setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
831              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
832    
833  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
834  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
835    
   
   
836  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),  setMethod("solve", signature(a = "diagonalMatrix", b = "missing"),
837            function(a, b, ...) {            function(a, b, ...) {
838                a@x <- 1/ a@x                a@x <- 1/ a@x
# Line 515  Line 841 
841            })            })
842    
843  solveDiag <- function(a, b, ...) {  solveDiag <- function(a, b, ...) {
844      if((n <- a@Dim[1]) != nrow(b))      if(a@Dim[1] != nrow(b))
845          stop("incompatible matrix dimensions")          stop("incompatible matrix dimensions")
846      ## trivially invert a 'in place' and multiply:      ## trivially invert a 'in place' and multiply:
847      a@x <- 1/ a@x      a@x <- 1/ a@x
# Line 534  Line 860 
860  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------  ###---------------- <Ops> (<Arith>, <Logic>, <Compare> ) ----------------------
861    
862  ## 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.)  
863  diagOdiag <- function(e1,e2) {  diagOdiag <- function(e1,e2) {
864      ## result should also be diagonal _ if possible _      ## result should also be diagonal _ if possible _
865      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 542  Line 867 
867      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,      r00 <- callGeneric(if(is.numeric(e1@x)) 0 else FALSE,
868                         if(is.numeric(e2@x)) 0 else FALSE)                         if(is.numeric(e2@x)) 0 else FALSE)
869      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* diagonal
870          if(is.numeric(r)) {          if(is.numeric(r)) { # "double" *or* "integer"
871              if(is.numeric(e2@x)) {              if(is.numeric(e2@x)) {
872                  e2@x <- r; return(.diag.2N(e2)) }                  e2@x <- r; return(.diag.2N(e2)) }
873              if(!is.numeric(e1@x))              if(!is.numeric(e1@x))
874                  ## e.g. e1, e2 are logical;                  ## e.g. e1, e2 are logical;
875                  e1 <- as(e1, "dMatrix")                  e1 <- as(e1, "dMatrix")
876                if(!is.double(r)) r <- as.double(r)
877          }          }
878          else if(is.logical(r))          else if(is.logical(r))
879              e1 <- as(e1, "lMatrix")              e1 <- as(e1, "lMatrix")
880          else stop("intermediate 'r' is of type", typeof(r))          else stop(gettextf("intermediate 'r' is of type %s",
881                               typeof(r)), domain=NA)
882          e1@x <- r          e1@x <- r
883          .diag.2N(e1)          .diag.2N(e1)
884      }      }
885      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
886            ## e.g., m == m
887          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
888          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
889            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         if(getOption("verbose"))  
             message("exploding  <diag>  o  <diag>  into dense matrix")  
890          d <- e1@Dim          d <- e1@Dim
891          n <- d[1]          n <- d[1]
892          stopifnot(length(r) == n)          stopifnot(length(r) == n)
893            if(isNum && !is.double(r)) r <- as.double(r)
894          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))          xx <- as.vector(matrix(rbind(r, matrix(r00,n,n)), n,n))
895          newcl <-          newcl <-
896              paste(if(isNum) "d" else if(isLog) {              paste0(if(isNum) "d" else if(isLog) {
897                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"                  if(!any(is.na(r)) && !any(is.na(r00))) "n" else "l"
898              } else stop("not yet implemented .. please report")              } else stop("not yet implemented .. please report"), "syMatrix")
                   ,  
                   "syMatrix", sep='')  
899    
900          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)          new(newcl, Dim = e1@Dim, Dimnames = e1@Dimnames, x = xx)
901      }      }
# Line 588  Line 913 
913              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
914  }  }
915    
916  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
917  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
918  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
919  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
920    diagOtri <- function(e1,e2) {
921        ## result must be triangular
922        r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
923        ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
924        e1.0 <- if(is.numeric(d1)) 0 else FALSE
925        r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
926        if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
927            diag(e2) <- r
928            ## check what happens "in the triangle"
929            e2.2 <- if(.n2) 2 else TRUE
930            if(!callGeneric(e1.0, e2.2) == e2.2) { # values "in triangle" can change:
931                n <- dim(e2)[1L]
932                it <- indTri(n, upper = (e2@uplo == "U"))
933                e2[it] <- callGeneric(e1.0, e2[it])
934            }
935            e2
936        }
937        else { ## result not triangular ---> general
938            rr <- as(e2, "generalMatrix")
939            diag(rr) <- r
940            rr
941        }
942    }
943    
944    
945    setMethod("Ops", signature(e1 = "diagonalMatrix", e2 = "triangularMatrix"),
946              diagOtri)
947    ## For the reverse,  Ops == "Arith" | "Compare" | "Logic"
948    ##   'Arith'  :=  '"+"', '"-"', '"*"', '"^"', '"%%"', '"%/%"', '"/"'
949    setMethod("Arith", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
950              function(e1,e2)
951          { ## this must only trigger for *dense* e1
952              switch(.Generic,
953                     "+" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"),   .diag.x(e2)),
954                     "-" = .Call(dtrMatrix_addDiag, as(e1,"dtrMatrix"), - .diag.x(e2)),
955                     "*" = {
956                         n <- e2@Dim[1L]
957                         d2 <- if(e2@diag == "U") { # unit-diagonal
958                             d <- rep.int(as1(e2@x), n)
959                             e2@x <- d
960                             e2@diag <- "N"
961                             d
962                         } else e2@x
963                         e2@x <- diag(e1) * d2
964                         e2
965                     },
966                     "^" = { ## will be dense ( as  <ANY> ^ 0 == 1 ):
967                         e1 ^ as(e2, "denseMatrix")
968                     },
969                     ## otherwise:
970                     callGeneric(e1, diag2Tsmart(e2,e1)))
971    })
972    
973    ## Compare --> 'swap' (e.g.   e1 < e2   <==>  e2 > e1 ):
974    setMethod("Compare", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
975              .Cmp.swap)
976    ## '&' and "|'  are commutative:
977    setMethod("Logic", signature(e1 = "triangularMatrix", e2 = "diagonalMatrix"),
978              function(e1,e2) callGeneric(e2, e1))
979    
980  ## For almost everything else, diag* shall be treated "as sparse" :  ## For almost everything else, diag* shall be treated "as sparse" :
981  ## These are cheap implementations via coercion  ## These are cheap implementations via coercion
# Line 615  Line 999 
999  ##       Compare --> logical  ##       Compare --> logical
1000  ##       Logic   --> logical: "lMatrix"  ##       Logic   --> logical: "lMatrix"
1001    
1002  ##  other = "numeric" : stay diagonal if possible  ## Other = "numeric" : stay diagonal if possible
1003  ## ddi*: Arith: result numeric, potentially ddiMatrix  ## ddi*: Arith: result numeric, potentially ddiMatrix
1004  setMethod("Arith", signature(e1 = "ddiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1005    setMethod("Arith", signature(e1 = "ddiMatrix", e2 = arg2),
1006            function(e1,e2) {            function(e1,e2) {
1007                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1008                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1009                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1010                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1011                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    if(e1@diag == "U") {
1012                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
1013                        e1@diag <- "N"                        e1@diag <- "N"
1014                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1015                    } else                        } ## else: result = e1  (is "U" diag)
1016                      } else {
1017                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1018                    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
1019                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1020                }                }
1021                      e1
1022                  } else
1023                callGeneric(diag2tT.u(e1,e2, "d"), e2)                callGeneric(diag2tT.u(e1,e2, "d"), e2)
1024            })            })
1025    
1026  setMethod("Arith", signature(e1 = "numeric", e2 = "ddiMatrix"),  for(arg1 in c("numeric","logical"))
1027    setMethod("Arith", signature(e1 = arg1, e2 = "ddiMatrix"),
1028            function(e1,e2) {            function(e1,e2) {
1029                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
1030                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1031                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1032                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1033                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    if(e2@diag == "U") {
1034                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1035                        e2@diag <- "N"                        e2@diag <- "N"
1036                        if(L1) r <- rep.int(r, n)                            e2@x[seq_len(n)] <- r # possibly recycling r
1037                    } else                        } ## else: result = e2  (is "U" diag)
1038                      } else {
1039                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1040                    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
1041                    return(e2)                        e2@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1042                }                }
1043                      e2
1044                  } else
1045                callGeneric(e1, diag2tT.u(e2,e1, "d"))                callGeneric(e1, diag2tT.u(e2,e1, "d"))
1046            })            })
1047    
1048  ## ldi* Arith --> result numeric, potentially ddiMatrix  ## ldi* Arith --> result numeric, potentially ddiMatrix
1049  setMethod("Arith", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1050    setMethod("Arith", signature(e1 = "ldiMatrix", e2 = arg2),
1051            function(e1,e2) {            function(e1,e2) {
1052                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1053                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1054                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1055                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1056                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1057                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "double"
1058                        e1@diag <- "N"                    if(e1@diag == "U") {
1059                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1060                    } else                            E@diag <- "N"
1061                              E@x[seq_len(n)] <- r # possibly recycling r
1062                          } ## else: result = E  (is "U" diag)
1063                      } else {
1064                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1065                    e1 <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1066                    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)  
1067                }                }
1068                callGeneric(diag2tT.u(e1,e2, "d"), e2)                    E
1069                  } else
1070                      callGeneric(diag2tT.u(e1,e2, "l"), e2)
1071            })            })
1072    
1073  setMethod("Arith", signature(e1 = "numeric", e2 = "ldiMatrix"),  for(arg1 in c("numeric","logical"))
1074    setMethod("Arith", signature(e1 = arg1, e2 = "ldiMatrix"),
1075            function(e1,e2) {            function(e1,e2) {
1076                n <- e2@Dim[1]; nsq <- n*n                n <- e2@Dim[1]
1077                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1078                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1079                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1080                    if(!L1 && le != nsq) e1 <- rep(e1, length.out = nsq)                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1081                    if(e2@diag == "U" && any((r <- callGeneric(e1, 1)) != 1)) {                    ## storage.mode(E@x) <- "double"
1082                        e2@diag <- "N"                    if(e2@diag == "U") {
1083                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(e1, 1)) != 1)) {
1084                    } else                            E@diag <- "N"
1085                              E@x[seq_len(n)] <- r # possibly recycling r
1086                          } ## else: result = E  (is "U" diag)
1087                      } else {
1088                        r <- callGeneric(e1, e2@x)                        r <- callGeneric(e1, e2@x)
1089                    e2 <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1090                    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)  
1091                }                }
1092                callGeneric(e1, diag2tT.u(e2,e1, "d"))                    E
1093                  } else
1094                      callGeneric(e1, diag2tT.u(e2,e1, "l"))
1095            })            })
1096    
1097  ## ddi*: for "Ops" without Arith --> result logical, potentially ldi  ## ddi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1098  setMethod("Ops", signature(e1 = "ddiMatrix", e2 = "numeric"),  ##
1099    ## Note that  ("numeric", "ddiMatrix")  is simply swapped, e.g.,
1100    if(FALSE) {
1101        selectMethod("<", c("numeric","lMatrix"))# Compare
1102        selectMethod("&", c("numeric","lMatrix"))# Logic
1103    } ## so we don't need to define a method here :
1104    
1105    for(arg2 in c("numeric","logical"))
1106    setMethod("Ops", signature(e1 = "ddiMatrix", e2 = arg2),
1107            function(e1,e2) {            function(e1,e2) {
1108                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1109                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1110                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1111                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1112                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
1113                    if(e1@diag == "U" && any((r <- callGeneric(1, e2)) != 1)) {                    ## storage.mode(E@x) <- "logical"
1114                        e1@diag <- "N"                    if(e1@diag == "U") {
1115                        if(L1) r <- rep.int(r, n)                        if(any((r <- callGeneric(1, e2)) != 1)) {
1116                    } else                            E@diag <- "N"
1117                              E@x[seq_len(n)] <- r # possibly recycling r
1118                          } ## else: result = E  (is "U" diag)
1119                      } else {
1120                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1121                    e1 <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))                        ## "future fixme": if we have idiMatrix, and r is 'integer', use idiMatrix
1122                    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)  
1123                }                }
1124                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    E
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ddiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               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)  
1125                    } else                    } else
1126                        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"))  
1127            })            })
1128    
1129  ## ldi*: for "Ops" without Arith --> result logical, potentially ldi  ## ldi*: for "Ops" without "Arith": <Compare> or <Logic> --> result logical, potentially ldi
1130  setMethod("Ops", signature(e1 = "ldiMatrix", e2 = "numeric"),  for(arg2 in c("numeric","logical"))
1131    setMethod("Ops", signature(e1 = "ldiMatrix", e2 = arg2),
1132            function(e1,e2) {            function(e1,e2) {
1133                n <- e1@Dim[1]; nsq <- n*n                n <- e1@Dim[1]
1134                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1135                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1136                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1137                    if(!L1 && le != nsq) e2 <- rep(e2, length.out = nsq)  
1138                    if(e1@diag == "U" && any((r <- callGeneric(TRUE, e2)) != 1)) {                    if(e1@diag == "U") {
1139                          if(any((r <- callGeneric(TRUE, e2)) != 1)) {
1140                        e1@diag <- "N"                        e1@diag <- "N"
1141                        if(L1) r <- rep.int(r, n)                            e1@x[seq_len(n)] <- r # possibly recycling r
1142                    } else                        } ## else: result = e1  (is "U" diag)
1143                      } else {
1144                        r <- callGeneric(e1@x, e2)                        r <- callGeneric(e1@x, e2)
1145                    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
1146                    return(e1)                        e1@x[] <- if(L1) r else r[1L + ((n+1)*(0:(n-1L))) %% le]
1147                }                }
1148                callGeneric(diag2tT.u(e1,e2, "l"), e2)                    e1
           })  
   
 setMethod("Ops", signature(e1 = "numeric", e2 = "ldiMatrix"),  
           function(e1,e2) {  
               n <- e2@Dim[1]; nsq <- n*n  
               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)  
1149                    } else                    } else
1150                        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"))  
1151            })            })
1152    
1153    
   
1154  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }  ## Not {"sparseMatrix", "numeric} :  {"denseMatrix", "matrix", ... }
1155  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
1156      ## ddi*:      ## ddi*:
1157      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
1158                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
1159      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
1160                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
1161      ## ldi*:      ## ldi*:
1162      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
1163                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
1164      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
1165                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
1166  }  }
1167  ## This should *not* dispatch to <dense> methods (in ./Ops.R ), as  
1168  ##  FALSE & <anything> |-> FALSE : hence result should be diagonal:  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1169  for(cl in diCls) {  if(FALSE) # now also contains "geMatrix"
1170      setMethod("&", signature(e1 = cl, e2 = "ANY"),  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1171                function(e1,e2) e1 & as(e2,"Matrix"))                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1172      setMethod("&", signature(e1 = "ANY", e2 = cl),  dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1173                function(e1,e2) as(e1,"Matrix") & e2)  for(DI in diCls) {
1174      for(c2 in c("denseMatrix", "Matrix")) {      dMeth <- if(extends(DI, "dMatrix"))
1175          setMethod("&", signature(e1 = cl, e2 = c2),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
1176                    function(e1,e2) e1 & Diagonal(x = diag(e2)))      else # "lMatrix", the only other kind for now
1177          setMethod("&", signature(e1 = c2, e2 = cl),          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2)
1178                    function(e1,e2) Diagonal(x = diag(e1)) & e2)      for(c2 in c(dense.subCl, "Matrix")) {
1179            for(Fun in c("*", "&")) {
1180                setMethod(Fun, signature(e1 = DI, e2 = c2),
1181                          function(e1,e2) callGeneric(e1, Diagonal(x = diag(e2))))
1182                setMethod(Fun, signature(e1 = c2, e2 = DI),
1183                          function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1184            }
1185            setMethod("^", signature(e1 = c2, e2 = DI),
1186                      function(e1,e2) callGeneric(Diagonal(x = diag(e1)), e2))
1187            for(Fun in c("%%", "%/%", "/")) ## 0 <op> 0 |--> NaN  for these.
1188                setMethod(Fun, signature(e1 = DI, e2 = c2), dMeth)
1189      }      }
1190  }  }
1191    
1192    ## Group methods "Math", "Math2" in                     --> ./Math.R
1193    
1194  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1195  ### ----------  any, all: separately here  ### ----------   the last 4: separately here
1196  for(cl in diCls) {  for(cl in diCls) {
1197  setMethod("any", cl,  setMethod("any", cl,
1198            function (x, ..., na.rm) {            function (x, ..., na.rm) {
1199                if(any(x@Dim == 0)) FALSE                if(any(x@Dim == 0)) FALSE
1200                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)
1201            })            })
1202  setMethod("all",  cl, function (x, ..., na.rm) any(x@Dim == 0))  setMethod("all",  cl, function (x, ..., na.rm) {
1203  setMethod("prod", cl, function (x, ..., na.rm) as.numeric(any(x@Dim == 0)))      n <- x@Dim[1]
1204        if(n >= 2) FALSE
1205        else if(n == 0 || x@diag == "U") TRUE
1206        else all(x@x, ..., na.rm = na.rm)
1207    })
1208    setMethod("prod", cl, function (x, ..., na.rm) {
1209        n <- x@Dim[1]
1210        if(n >= 2) 0
1211        else if(n == 0 || x@diag == "U") 1
1212        else ## n == 1, diag = "N" :
1213            prod(x@x, ..., na.rm = na.rm)
1214    })
1215    
1216  setMethod("sum", cl,  setMethod("sum", cl,
1217            function(x, ..., na.rm) {            function(x, ..., na.rm) {
# Line 847  Line 1250 
1250      invisible(x)      invisible(x)
1251  }  }
1252    
1253    ## somewhat consistent with "print" for sparseMatrix :
1254    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
1255    
1256  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
1257            function(object) {            function(object) {
1258                d <- dim(object)                d <- dim(object)
1259                cl <- class(object)                cl <- class(object)
1260                cat(sprintf('%d x %d diagonal matrix of class "%s"\n',                cat(sprintf('%d x %d diagonal matrix of class "%s"',
1261                            d[1], d[2], cl))                            d[1], d[2], cl))
1262                  if(d[1] < 50) {
1263                      cat("\n")
1264                prDiag(object)                prDiag(object)
1265                  } else {
1266                      cat(", with diagonal entries\n")
1267                      show(diag(object))
1268                      invisible(object)
1269                  }
1270              })
1271    
1272    rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1273       dense.subCl, diCls)# not used elsewhere
1274    
1275    setMethod("summary", signature(object = "diagonalMatrix"),
1276              function(object, ...) {
1277                  d <- dim(object)
1278                  r <- summary(object@x, ...)
1279                  attr(r, "header") <-
1280                      sprintf('%d x %d diagonal Matrix of class "%s"',
1281                              d[1], d[2], class(object))
1282                  ## use ole' S3 technology for such a simple case
1283                  class(r) <- c("diagSummary", class(r))
1284                  r
1285            })            })
1286    
1287    print.diagSummary <- function (x, ...) {
1288        cat(attr(x, "header"),"\n")
1289        class(x) <- class(x)[-1]
1290        print(x, ...)
1291        invisible(x)
1292    }

Legend:
Removed from v.2197  
changed lines
  Added in v.3069

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