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

revision 2877, Wed May 8 20:59:04 2013 UTC revision 3079, Tue Mar 31 15:29:43 2015 UTC
# Line 120  Line 120 
120      ## block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
121      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
122    
123      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"      Tlst <- lapply(lapply(lst, as_Csp2), # includes "diagU2N"
124                     as, "TsparseMatrix")                     as, "TsparseMatrix")
125      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
126      ## else      ## else
# Line 198  Line 198 
198      as(.bdiag(alis), "CsparseMatrix")      as(.bdiag(alis), "CsparseMatrix")
199  }  }
200    
201    setMethod("tril", "diagonalMatrix", function(x, k = 0, ...)
202        if(k >= 0) x else .setZero(x, paste0(.M.kind(x), "tCMatrix")))
203    
204    setMethod("triu", "diagonalMatrix", function(x, k = 0, ...)
205        if(k <= 0) x else  .setZero(x, paste0(.M.kind(x), "tCMatrix")))
206    
207    
208    
209  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {  .diag2tT <- function(from, uplo = "U", kind = .M.kind(from)) {
210      ## to triangular Tsparse      ## to triangular Tsparse
# Line 227  Line 234 
234                  n))                  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 251  Line 258 
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  ##_no_longer_ 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", diag2tT)  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  ##_no_longer_ 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", diag2tT)  setAs("ldiMatrix", "lsparseMatrix", di2tT)
278  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
279        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
280  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix", function(from) .diag2sT(from, "U", "l"))
281        function(from) .diag2sT(from, "U", "l"))  rm(di2tT)
   
282    
283  setAs("diagonalMatrix", "nMatrix",  setAs("diagonalMatrix", "nMatrix",
284        function(from) {        function(from) {
           n <- from@Dim[1]  
285            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L            i <- if(from@diag == "U") integer(0) else which(isN0(from@x)) - 1L
286            new("ntTMatrix", i = i, j = i, diag = from@diag,            new("ntTMatrix", i = i, j = i, diag = from@diag,
287                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
# Line 285  Line 289 
289    
290  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))  setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
291    
292  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ##' A version of diag(x,n) which *does* preserve the mode of x, where diag() "fails"
293  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
294      y <- matrix(as0(mod=mode(x)), n,n)      y <- matrix(as0(mod=mode(x)), n,n)
295      if (n > 0) y[1L + 0:(n - 1L) * (n + 1)] <- 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 515  Line 525 
525  setMethod("t", signature(x = "diagonalMatrix"),  setMethod("t", signature(x = "diagonalMatrix"),
526            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })            function(x) { x@Dimnames <- x@Dimnames[2:1] ; x })
527    
528  setMethod("isDiagonal", signature(object = "diagonalMatrix"),  setMethod("isDiagonal",   "diagonalMatrix", function(object) TRUE)
529            function(object) TRUE)  setMethod("isTriangular", "diagonalMatrix", function(object, upper=NA, ...) TRUE)
530  setMethod("isTriangular", signature(object = "diagonalMatrix"),  setMethod("isSymmetric",  "diagonalMatrix", function(object, ...) TRUE)
           function(object) TRUE)  
 setMethod("isSymmetric", signature(object = "diagonalMatrix"),  
           function(object, ...) TRUE)  
531    
532  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)  setMethod("symmpart", signature(x = "diagonalMatrix"), function(x) x)
533  setMethod("skewpart", signature(x = "diagonalMatrix"), setZero)  setMethod("skewpart", signature(x = "diagonalMatrix"), function(x) .setZero(x))
534    
535  setMethod("chol", signature(x = "ddiMatrix"),  setMethod("chol", signature(x = "ddiMatrix"),
536            function(x, pivot, ...) {            function(x, pivot, ...) {
# Line 570  Line 577 
577      return(y)      return(y)
578  }  }
579    
580    ##' Boolean Algebra/Arithmetic Product of Diagonal Matrices
581    ##'  %&%
582    diagdiagprodBool <- function(x, y) {
583        dimCheck(x,y)
584        if(x@diag != "U") {
585            if(!is.logical(x@x)) x <- as(x, "lMatrix")
586            if(y@diag != "U") {
587                nx <- x@x & y@x
588                x@x <- as.logical(nx)
589            }
590            x
591        } else { ## x is unit diagonal: return y
592            if(!is.logical(y@x)) y <- as(y, "lMatrix")
593            y
594        }
595    }
596    
597  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
598            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprod, valueClass = "ddiMatrix")
599    
600  formals(diagdiagprod) <- alist(x=, y=x)  setMethod("%&%", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
601              diagdiagprodBool, valueClass = "ldiMatrix")# do *not* have "ndiMatrix" !
602    
603    ##' Both Numeric or Boolean Algebra/Arithmetic Product of Diagonal Matrices
604    diagdiagprodFlexi <- function(x, y=NULL, boolArith = NA, ...)
605    {
606        dimCheck(x,y)
607        bool <- isTRUE(boolArith)
608        if(x@diag != "U") {
609            if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
610            if(y@diag != "U") {
611                if(bool) {
612                    nx <- x@x & y@x
613                    x@x <- as.logical(nx)
614                } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
615                    nx <- x@x * y@x
616                    if(is.numeric(nx) && !is.numeric(x@x))
617                        x <- as(x, "dMatrix")
618                    x@x <- as.numeric(nx)
619                }
620            }
621            x
622        } else { ## x is unit diagonal: return y
623            if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
624            y
625        }
626    }
627  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
628            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
629  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "diagonalMatrix"),
630            diagdiagprod, valueClass = "ddiMatrix")            diagdiagprodFlexi)
631  setMethod("crossprod", signature(x = "diagonalMatrix", y = "missing"),  
632            diagdiagprod, valueClass = "ddiMatrix")  ##' crossprod(x) := x'x
633  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"),  diagprod <- function(x, y = NULL, boolArith = NA, ...) {
634            diagdiagprod, valueClass = "ddiMatrix")      bool <- isTRUE(boolArith)
635        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
636        if(x@diag != "U") {
637            if(bool) {
638                nx <- x@x & y@x
639                x@x <- as.logical(nx)
640            } else { ## boolArith is NA or FALSE: ==> numeric, as have *no* "diagMatrix" patter[n]:
641                nx <- x@x * x@x
642                if(is.numeric(nx) && !is.numeric(x@x))
643                    x <- as(x, "dMatrix")
644                x@x <- as.numeric(nx)
645            }
646        }
647        x
648    }
649    setMethod( "crossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
650    setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "missing"), diagprod)
651    
652    
653    ## analogous to matdiagprod() below:
654  diagmatprod <- function(x, y) {  diagmatprod <- function(x, y) {
655      ## x is diagonalMatrix      ## x is diagonalMatrix
656      dx <- dim(x)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
657      dy <- dim(y)      Matrix(if(x@diag == "U") y else x@x * y)
658      if(dx[2] != dy[1]) stop("non-matching dimensions")  }
659      as(if(x@diag == "U") y else x@x * y, "Matrix")  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"), diagmatprod)
660  }  
661  setMethod("%*%", signature(x = "diagonalMatrix", y = "matrix"),  formals(diagmatprod) <- alist(x=, y=NULL, boolArith = NA, ...=) ## FIXME boolArith
662            diagmatprod)  diagmatprod2 <- function(x, y=NULL, boolArith = NA, ...) {
663  ## sneaky .. :      ## x is diagonalMatrix
664  formals(diagmatprod) <- alist(x=, y=NULL)      if(x@Dim[2] != nrow(y)) stop("non-matching dimensions")
665  setMethod("crossprod", signature(x = "diagonalMatrix", y = "matrix"),      Matrix(if(x@diag == "U") y else x@x * y)
666            diagmatprod)  }
667    setMethod("crossprod",  signature(x = "diagonalMatrix", y = "matrix"), diagmatprod2)
668    
669  diagGeprod <- function(x, y) {  diagGeprod <- function(x, y) {
670      dx <- dim(x)      if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
     dy <- dim(y)  
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
671      if(x@diag != "U")      if(x@diag != "U")
672          y@x <- x@x * y@x          y@x <- x@x * y@x
673      y      y
674  }  }
675  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "dgeMatrix"), diagGeprod)
676  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)  setMethod("%*%", signature(x= "diagonalMatrix", y= "lgeMatrix"), diagGeprod)
 formals(diagGeprod) <- alist(x=, y=NULL)  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"),  
           diagGeprod, valueClass = "dgeMatrix")  
 setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"),  
           diagGeprod)  
677    
678    diagGeprodBool <- function(x, y) {
679        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
680        if(!is.logical(y@x)) y <- as(y, "lMatrix")
681        if(x@diag != "U")
682            y@x <- x@x & y@x
683        y
684    }
685    setMethod("%&%", signature(x= "diagonalMatrix", y= "geMatrix"), diagGeprodBool)
686    
687    diagGeprod2 <- function(x, y=NULL, boolArith = NA, ...) {
688        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
689        bool <- isTRUE(boolArith)
690        if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
691        if(x@diag != "U")
692            y@x <- if(bool) x@x & y@x else x@x * y@x
693        y
694    }
695    setMethod("crossprod", signature(x = "diagonalMatrix", y = "dgeMatrix"), diagGeprod2)
696    setMethod("crossprod", signature(x = "diagonalMatrix", y = "lgeMatrix"), diagGeprod2)
697    
698    
699    ## analogous to diagmatprod() above:
700  matdiagprod <- function(x, y) {  matdiagprod <- function(x, y) {
701      dx <- dim(x)      dx <- dim(x)
702      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
703      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))      Matrix(if(y@diag == "U") x else x * rep(y@x, each = dx[1]))
704  }  }
705  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "matrix", y = "diagonalMatrix"), matdiagprod)
           matdiagprod)  
 formals(matdiagprod) <- alist(x=, y=NULL)  
 setMethod("tcrossprod", signature(x = "matrix", y = "diagonalMatrix"),  
           matdiagprod)  
706    
707  gediagprod <- function(x, y) {  gediagprod <- function(x, y) {
708      dx <- dim(x)      dx <- dim(x)
709      dy <- dim(y)      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
     if(dx[2] != dy[1]) stop("non-matching dimensions")  
710      if(y@diag == "N")      if(y@diag == "N")
711          x@x <- x@x * rep(y@x, each = dx[1])          x@x <- x@x * rep(y@x, each = dx[1])
712      x      x
713  }  }
714  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "dgeMatrix", y= "diagonalMatrix"), gediagprod)
715  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)  setMethod("%*%", signature(x= "lgeMatrix", y= "diagonalMatrix"), gediagprod)
716  formals(gediagprod) <- alist(x=, y=NULL)  
717  setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"),  gediagprodBool <- function(x, y) {
718            gediagprod)      dx <- dim(x)
719  setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"),      if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
720            gediagprod)      if(!is.logical(x@x)) x <- as(x, "lMatrix")
721        if(y@diag == "N")
722            x@x <- x@x & rep(y@x, each = dx[1])
723        x
724    }
725    setMethod("%&%", signature(x= "geMatrix", y= "diagonalMatrix"), gediagprodBool)
726    
727    setMethod("tcrossprod",signature(x = "matrix", y = "diagonalMatrix"),
728              function(x, y=NULL, boolArith = NA, ...) {
729                  dx <- dim(x)
730                  if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
731                  bool <- isTRUE(boolArith)
732                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
733                  Matrix(if(y@diag == "U") x else
734                         if(bool) x & rep(y@x, each = dx[1])
735                         else     x * rep(y@x, each = dx[1]))
736              })
737    
738    setMethod("crossprod", signature(x = "matrix", y = "diagonalMatrix"),
739              function(x, y=NULL, boolArith = NA, ...) {
740                  dx <- dim(x)
741                  if(dx[1] != y@Dim[1]) stop("non-matching dimensions")
742                  bool <- isTRUE(boolArith)
743                  if(bool && !is.logical(y@x)) y <- as(y, "lMatrix")
744                  Matrix(if(bool) t(rep.int(y@x, dx[2]) & x)
745                             else t(rep.int(y@x, dx[2]) * x))
746              })
747    
748    
749    gediagprod2 <- function(x, y=NULL, boolArith = NA, ...) {
750        dx <- dim(x)
751        if(dx[2] != y@Dim[1]) stop("non-matching dimensions")
752        bool <- isTRUE(boolArith)
753        if(bool && !is.logical(x@x)) x <- as(x, "lMatrix")
754        if(y@diag == "N")
755            x@x <- if(bool) x@x & rep(y@x, each = dx[1])
756                   else     x@x * rep(y@x, each = dx[1])
757        x
758    }
759    setMethod("tcrossprod", signature(x = "dgeMatrix", y = "diagonalMatrix"), gediagprod2)
760    setMethod("tcrossprod", signature(x = "lgeMatrix", y = "diagonalMatrix"), gediagprod2)
761    
762    
763  ## crossprod {more of these}  ## crossprod {more of these}
764    
# Line 657  Line 775 
775  ##        function(x, y = NULL) {  ##        function(x, y = NULL) {
776  ##           })  ##           })
777    
778  Cspdiagprod <- function(x, y) {  ##' @param x CsparseMatrix
779      dx <- dim(x <- .Call(Csparse_diagU2N, x))  ##' @param y diagonalMatrix
780      dy <- dim(y)  ##' @return x %*% y
781      if(dx[2] != dy[1]) stop("non-matching dimensions")  Cspdiagprod <- function(x, y, boolArith = NA, ...) {
782      if(y@diag == "N") {      if((m <- ncol(x)) != y@Dim[1]) stop("non-matching dimensions")
783          if(!all(y@x[1L] == y@x[-1L]) && is(x, "symmetricMatrix"))      if(y@diag == "N") { ## otherwise: y == Diagonal(n) : multiplication is identity
784            x <- .Call(Csparse_diagU2N, x)
785            cx <- getClass(class(x))
786            if(!all(y@x[1L] == y@x[-1L]) && extends(cx, "symmetricMatrix"))
787              x <- as(x, "generalMatrix")              x <- as(x, "generalMatrix")
788          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])
789            if(isTRUE(boolArith)) {
790                if(extends(cx, "nMatrix")) x <- as(x, "lMatrix") # so, has y@x
791                x@x <- r <- x@x & y@x[x@i + 1L]
792                if(!anyNA(r) && !extends(cx, "diagonalMatrix")) x <- as(drop0(x), "nMatrix")
793            } else {
794                if(!extends(cx, "dMatrix")) x <- as(x, "dMatrix") # <- FIXME if we have zMatrix
795          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
796      }      }
797      if(is(x, "compMatrix") && length(xf <- x@factors)) {          if(.hasSlot(x, "factors") && length(x@factors)) {# drop cashed ones
798          ## instead of dropping all factors, be smart about some          ## instead of dropping all factors, be smart about some
799          ## TODO ......          ## TODO ......
800          x@factors <- list()          x@factors <- list()
801      }      }
802      x      x
803        } else { #  y is unit-diagonal ==> "return x"
804            cx <- getClass(class(x))
805            if(isTRUE(boolArith)) {
806                is.l <- if(extends(cx, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
807                    x <- as(x, "lMatrix"); TRUE } else extends(cx, "lMatrix")
808                if(is.l && !anyNA(x@x)) as(drop0(x), "nMatrix")
809                else if(is.l) x else # defensive:
810                as(x, "lMatrix")
811            } else {
812                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
813                if(extends(cx, "nMatrix") || extends(cx, "lMatrix"))
814                    as(x, "dMatrix") else x
815            }
816        }
817  }  }
818    
819  diagCspprod <- function(x, y) {  ##' @param x diagonalMatrix
820      dx <- dim(x)  ##' @param y CsparseMatrix
821      dy <- dim(y <- .Call(Csparse_diagU2N, y))  ##' @return x %*% y
822      if(dx[2] != dy[1]) stop("non-matching dimensions")  diagCspprod <- function(x, y, boolArith = NA, ...) {
823        if(x@Dim[2] != y@Dim[1]) stop("non-matching dimensions")
824      if(x@diag == "N") {      if(x@diag == "N") {
825          if(!all(x@x[1L] == x@x[-1L]) && is(y, "symmetricMatrix"))          y <- .Call(Csparse_diagU2N, y)
826            cy <- getClass(class(y))
827            if(!all(x@x[1L] == x@x[-1L]) && extends(cy, "symmetricMatrix"))
828              y <- as(y, "generalMatrix")              y <- as(y, "generalMatrix")
829            if(isTRUE(boolArith)) {
830                if(extends(cy, "nMatrix")) y <- as(y, "lMatrix") # so, has y@x
831                y@x <- r <- y@x & x@x[y@i + 1L]
832                if(!anyNA(r) && !extends(cy, "diagonalMatrix")) y <- as(drop0(y), "nMatrix")
833            } else {
834                if(!extends(cy, "dMatrix")) y <- as(y, "dMatrix") # <- FIXME if we have zMatrix
835          y@x <- y@x * x@x[y@i + 1L]          y@x <- y@x * x@x[y@i + 1L]
836      }      }
837      if(is(y, "compMatrix") && length(yf <- y@factors)) {          if(.hasSlot(y, "factors") && length(y@factors)) {
838         ## if(.hasSlot(y, "factors") && length(yf <- y@factors)) { ## -- TODO? --
839          ## instead of dropping all factors, be smart about some          ## instead of dropping all factors, be smart about some
840          ## TODO              ## keep <- character()
841          keep <- character()              ## if(any(names(yf) == "LU")) { ## <- not easy: y = P'LUQ,  x y = xP'LUQ => LU ???
842          if(iLU <- names(yf) == "LU") {              ##     keep <- "LU"
843              ## TODO keep <- "LU"              ## }
844          }              ## y@factors <- yf[keep]
845          y@factors <- yf[keep]              y@factors <- list()
846      }      }
847      y      y
848        } else { ## x @ diag  == "U"
849            cy <- getClass(class(y))
850            if(isTRUE(boolArith)) {
851                is.l <- if(extends(cy, "dMatrix")) { ## <- FIXME: extend once we have iMatrix, zMatrix
852                    y <- as(y, "lMatrix"); TRUE } else extends(cy, "lMatrix")
853                if(is.l && !anyNA(y@x)) as(drop0(y), "nMatrix")
854                else if(is.l) y else # defensive:
855                as(y, "lMatrix")
856            } else {
857                ## else boolArith is  NA or FALSE {which are equivalent here, das diagonal = "numLike"}
858                if(extends(cy, "nMatrix") || extends(cy, "lMatrix"))
859                    as(y, "dMatrix") else y
860            }
861        }
862  }  }
863    
864    ## + 'boolArith' argument  { ==> .local() is used in any case; keep formals simple :}
865  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
866            function(x, y = NULL) diagCspprod(x, y))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, y, boolArith=boolArith))
867    
868  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("crossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
869            function(x, y = NULL) diagCspprod(x, as(y, "CsparseMatrix")))            function(x, y=NULL, boolArith=NA, ...)
870                  diagCspprod(x, as(y, "CsparseMatrix"), boolArith=boolArith))
871    
872  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway  ## Prefer calling diagCspprod to Cspdiagprod if going to transpose anyway
873  ##  x'y == (y'x)'  ##  x'y == (y'x)'
874  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
875            function(x, y = NULL) t(diagCspprod(y, x)))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, x, boolArith=boolArith)))
876    
877  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("crossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
878            function(x, y = NULL) t(diagCspprod(y, as(x, "Csparsematrix"))))            function(x, y=NULL, boolArith=NA, ...) t(diagCspprod(y, as(x, "Csparsematrix"), boolArith=boolArith)))
879    
880  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
881            function(x, y = NULL) diagCspprod(x, t(y)))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(y), boolArith=boolArith))
882    
883  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("tcrossprod", signature(x = "diagonalMatrix", y = "sparseMatrix"),
884            function(x, y = NULL) diagCspprod(x, t(as(y, "CsparseMatrix"))))            function(x, y=NULL, boolArith=NA, ...) diagCspprod(x, t(as(y, "CsparseMatrix")), boolArith=boolArith))
885    
886  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
887            function(x, y = NULL) Cspdiagprod(x, y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(x, y, boolArith=boolArith))
888    
889  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("tcrossprod", signature(x = "sparseMatrix", y = "diagonalMatrix"),
890            function(x, y = NULL) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y=NULL, boolArith=NA, ...) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=boolArith))
891    
892  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
893            function(x, y) diagCspprod(x, y))            function(x, y) diagCspprod(x, y, boolArith=NA))
894    setMethod("%&%", signature(x = "diagonalMatrix", y = "CsparseMatrix"),
895              function(x, y) diagCspprod(x, y, boolArith=TRUE))
896    
897    ## instead of "sparseMatrix", use: [RT]sparse.. ("closer" in method dispatch)
898    for(cl in c("TsparseMatrix", "RsparseMatrix")) {
899    
900  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),  setMethod("%*%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
901            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y))            function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=NA))
902    
903  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
904            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y))            function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=NA))
905    
906    setMethod("%&%", signature(x = "diagonalMatrix", y = "sparseMatrix"),
907              function(x, y) diagCspprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
908    
909    setMethod("%&%", signature(x = "sparseMatrix", y = "diagonalMatrix"),
910              function(x, y) Cspdiagprod(as(x, "CsparseMatrix"), y, boolArith=TRUE))
911    
912    }
913    
914  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),  setMethod("%*%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
915            function(x, y) Cspdiagprod(x, y))            function(x, y) Cspdiagprod(x, y, boolArith=NA))
916    setMethod("%&%", signature(x = "CsparseMatrix", y = "diagonalMatrix"),
917              function(x, y) Cspdiagprod(x, y, boolArith=TRUE))
918    
919  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*  ## TODO: Write tests in ./tests/ which ensure that many "ops" with diagonal*
920  ##       do indeed work by going through sparse (and *not* ddense)!  ##       do indeed work by going through sparse (and *not* ddense)!
# Line 825  Line 1007 
1007      ## result must be triangular      ## result must be triangular
1008      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"      r <- callGeneric(d1 <- .diag.x(e1), diag(e2)) # error if not "compatible"
1009      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:      ## Check what happens with non-diagonals, i.e. (0 o 0), (FALSE o 0), ...:
1010      e1.0 <- if(.n1 <- is.numeric(d1   )) 0 else FALSE      e1.0 <- if(is.numeric(d1)) 0 else FALSE
1011      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)      r00 <- callGeneric(e1.0, if(.n2 <- is.numeric(e2[0])) 0 else FALSE)
1012      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular      if(is0(r00)) { ##  r00 == 0 or FALSE --- result *is* triangular
1013          diag(e2) <- r          diag(e2) <- r
# Line 910  Line 1092 
1092            function(e1,e2) {            function(e1,e2) {
1093                n <- e1@Dim[1]                n <- e1@Dim[1]
1094                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1095                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1096                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1097                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1098                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 932  Line 1114 
1114            function(e1,e2) {            function(e1,e2) {
1115                n <- e2@Dim[1]                n <- e2@Dim[1]
1116                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1117                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1118                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1119                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1120                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 955  Line 1137 
1137            function(e1,e2) {            function(e1,e2) {
1138                n <- e1@Dim[1]                n <- e1@Dim[1]
1139                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1140                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1141                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1142                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ddiMatrix", check=FALSE)  
1143                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1144                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1145                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 980  Line 1161 
1161            function(e1,e2) {            function(e1,e2) {
1162                n <- e2@Dim[1]                n <- e2@Dim[1]
1163                f0 <- callGeneric(e1, 0)                f0 <- callGeneric(e1, 0)
1164                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1165                    L1 <- (le <- length(e1)) == 1L                    L1 <- (le <- length(e1)) == 1L
1166                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e2, "ddiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e2, "ddiMatrix", check=FALSE)  
1167                    ## storage.mode(E@x) <- "double"                    ## storage.mode(E@x) <- "double"
1168                    if(e2@diag == "U") {                    if(e2@diag == "U") {
1169                        if(any((r <- callGeneric(e1, 1)) != 1)) {                        if(any((r <- callGeneric(e1, 1)) != 1)) {
# Line 1013  Line 1193 
1193            function(e1,e2) {            function(e1,e2) {
1194                n <- e1@Dim[1]                n <- e1@Dim[1]
1195                f0 <- callGeneric(0, e2)                f0 <- callGeneric(0, e2)
1196                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1197                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1198                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"))#FIXME: if ok, check=FALSE                    E <- copyClass(e1, "ldiMatrix", c("diag", "Dim", "Dimnames"), check=FALSE)
                   ## E <- copyClass(e1, "ldiMatrix", check=FALSE)  
1199                    ## storage.mode(E@x) <- "logical"                    ## storage.mode(E@x) <- "logical"
1200                    if(e1@diag == "U") {                    if(e1@diag == "U") {
1201                        if(any((r <- callGeneric(1, e2)) != 1)) {                        if(any((r <- callGeneric(1, e2)) != 1)) {
# Line 1039  Line 1218 
1218            function(e1,e2) {            function(e1,e2) {
1219                n <- e1@Dim[1]                n <- e1@Dim[1]
1220                f0 <- callGeneric(FALSE, e2)                f0 <- callGeneric(FALSE, e2)
1221                if(all(is0(f0))) { # remain diagonal                if(all0(f0)) { # remain diagonal
1222                    L1 <- (le <- length(e2)) == 1L                    L1 <- (le <- length(e2)) == 1L
1223    
1224                    if(e1@diag == "U") {                    if(e1@diag == "U") {
# Line 1073  Line 1252 
1252  }  }
1253    
1254  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
1255    if(FALSE) # now also contains "geMatrix"
1256  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
1257                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })                         names(dM.scl)[vapply(dM.scl, slot, 0, "distance") == 1] })
1258    dense.subCl <- paste0(c("d","l","n"), "denseMatrix")
1259  for(DI in diCls) {  for(DI in diCls) {
1260      dMeth <- if(extends(DI, "dMatrix"))      dMeth <- if(extends(DI, "dMatrix"))
1261          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)          function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2)
# Line 1094  Line 1275 
1275      }      }
1276  }  }
1277    
1278    ## Group methods "Math", "Math2" in                     --> ./Math.R
1279    
1280  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"  ### "Summary" : "max"   "min"   "range" "prod"  "sum"   "any"   "all"
1281  ### ----------   the last 4: separately here  ### ----------   the last 4: separately here
# Line 1173  Line 1355 
1355                }                }
1356            })            })
1357    
1358  rm(dense.subCl, diCls)# not used elsewhere  rm(arg1, arg2, other, DI, Fun, cl, c1, c2,
1359       dense.subCl, diCls)# not used elsewhere
1360    
1361  setMethod("summary", signature(object = "diagonalMatrix"),  setMethod("summary", signature(object = "diagonalMatrix"),
1362            function(object, ...) {            function(object, ...) {

Legend:
Removed from v.2877  
changed lines
  Added in v.3079

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge