# SCM Repository

[matrix] Diff of /pkg/Matrix/R/diagMatrix.R
 [matrix] / pkg / Matrix / R / diagMatrix.R

# Diff of /pkg/Matrix/R/diagMatrix.R

revision 2984, Sat Apr 12 21:37:37 2014 UTC revision 2998, Tue Aug 12 10:27:15 2014 UTC
# Line 227  Line 227
227                  n))                  n))
228  }  }
229
230  ## diagonal -> triangular,  upper / lower depending on "partner":  ## diagonal -> triangular,  upper / lower depending on "partner" 'x':
231  diag2tT.u <- function(d, x, kind = .M.kind(d))  diag2tT.u <- function(d, x, kind = .M.kind(d))
232      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)      .diag2tT(d, uplo = if(is(x,"triangularMatrix")) x@uplo else "U", kind)
233
# Line 259  Line 259
259  setAs("ddiMatrix", "dsparseMatrix", di2tT)  setAs("ddiMatrix", "dsparseMatrix", di2tT)
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", function(from) .diag2sT(from, "U", "d"))
function(from) .diag2sT(from, "U", "d"))
263  ##  ##
264  ## ldi*:  ## ldi*:
265  di2tT <- function(from) .diag2tT(from, "U", "l")  di2tT <- function(from) .diag2tT(from, "U", "l")
# Line 271  Line 270
270  setAs("ldiMatrix", "lsparseMatrix", di2tT)  setAs("ldiMatrix", "lsparseMatrix", di2tT)
271  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
272        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
273  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
function(from) .diag2sT(from, "U", "l"))
274  rm(di2tT)  rm(di2tT)
275
276  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
# Line 285  Line 283
283
284  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
285
286  ## 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"
287  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
288      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
289      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- x
290      y      y
291  }  }
292    ## NB: diag(x,n) is really faster for n >= 20, and even more for large n
293    ## --> using diag() where possible, ==> .ddi2mat()
294
295  setAs("diagonalMatrix", "matrix",  .diag2mat <- function(from)
296        function(from) {      ## want "ldiMatrix" -> <logical> "matrix"  (but integer -> <double> for now)
297            ## want "ldiMatrix" -> <logical> "matrix" :      mkDiag(if(from@diag == "U") as1(from@x) else from@x, n = from@Dim[1])
298            mkDiag(if(from@diag == "U") as1(from@x) else from@x,
299                   n = from@Dim[1])  .ddi2mat <- function(from)
300        })      base::diag(if(from@diag == "U") as1(from@x) else from@x, nrow = from@Dim[1])
301
302    setAs("ddiMatrix", "matrix", .ddi2mat)
303    ## the non-ddi diagonalMatrix -- only "ldiMatrix" currently:
304    setAs("diagonalMatrix", "matrix", .diag2mat)
305
306  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),  setMethod("as.vector", signature(x = "diagonalMatrix", mode="missing"),
307            function(x, mode) {            function(x, mode) {
# Line 581  Line 585
585            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
586
587
588    ## analogous to matdiagprod() below:
589  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
590      ## x is diagonalMatrix      ## x is diagonalMatrix
591      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
592      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
if(dx[2] != dy[1]) stop("non-matching dimensions")
as(if(x@diag == "U") y else x@x * y, "Matrix")
593  }  }
594  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),
595            diagmatprod)            diagmatprod)
596  ## sneaky .. :  ## sneaky .. :
597  formals(diagmatprod) <- alist(x=, y=NULL)  formals(diagmatprod) <- alist(x=, y=NULL)
598  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),  setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
diagmatprod)
599
600  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
601      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")
602      if(x@diag != "U")      if(x@diag != "U")
603          y@x <- x@x * y@x          y@x <- x@x * y@x
604      y      y
# Line 611  Line 611
611  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),
612            diagGeprod)            diagGeprod)
613
614    ## analogous to diagmatprod() above:
615  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
616      dx <- dim(x)      dx <- dim(x)
617      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
if(dx[2] != dy[1]) stop("non-matching dimensions")
618      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]))
619  }  }
620  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
matdiagprod)
621  formals(matdiagprod) <- alist(x=, y=NULL)  formals(matdiagprod) <- alist(x=, y=NULL)
622  setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
623            matdiagprod)  setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
624              function(x, y=NULL) {
625                  dx <- dim(x)
626                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
627                  Matrix(t(rep.int(y@x, dx[2]) * x))
628              })
629
630  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
631      dx <- dim(x)      dx <- dim(x)
632      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
if(dx[2] != dy[1]) stop("non-matching dimensions")
633      if(y@diag == "N")      if(y@diag == "N")
634          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
635      x      x
# Line 658  Line 661
661  ##' @param y diagonalMatrix  ##' @param y diagonalMatrix
662  ##' @return x %*% y  ##' @return x %*% y
663  Cspdiagprod <- function(x, y) {  Cspdiagprod <- function(x, y) {
664      dx <- dim(x <- .Call(Csparse_diagU2N, x))      m <- ncol(x <- .Call(Csparse_diagU2N, x))
665      dy <- dim(y)      if(m != y@Dim[1]) stop("non-matching dimensions")
if(dx[2] != dy[1]) stop("non-matching dimensions")
666      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
667          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))
668              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
669          ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])          ind <- rep.int(seq_len(m), x@p[-1] - x@p[-m-1L])
670          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
671          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
# Line 679  Line 681
681  ##' @param y CsparseMatrix  ##' @param y CsparseMatrix
682  ##' @return x %*% y  ##' @return x %*% y
683  diagCspprod <- function(x, y) {  diagCspprod <- function(x, y) {
684      dx <- dim(x)      y <- .Call(Csparse_diagU2N, y)
685      dy <- dim(y <- .Call(Csparse_diagU2N, y))      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
if(dx[2] != dy[1]) stop("non-matching dimensions")
686      if(x@diag == "N") {      if(x@diag == "N") {
687          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))
688              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
# Line 1099  Line 1100
1100      }      }
1101  }  }
1102

1103  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1104  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
1105  for(cl in diCls) {  for(cl in diCls) {
# Line 1178  Line 1178
1178                }                }
1179            })            })
1180
1181  rm(arg1, arg2, other, DI, cl, c1, c2,  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1182     dense.subCl, diCls)# not used elsewhere     dense.subCl, diCls)# not used elsewhere
1183
1184  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),

Legend:
 Removed from v.2984 changed lines Added in v.2998