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 2458, Fri Aug 28 22:19:48 2009 UTC pkg/Matrix/R/diagMatrix.R revision 2608, Tue Aug 10 08:23:48 2010 UTC
# Line 32  Line 32 
32      }      }
33  }  }
34    
35  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {  .sparseDiagonal <- function(n, x = rep.int(1,m), uplo = "U",
36                                shape = if(missing(cols)) "t" else "g",
37                                kind, cols = if(n) 0:(n - 1L) else integer(0))
38    {
39      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
40      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if(!missing(cols))
41      else if(lx != n) stop("length(x) must be 1 or n")          stopifnot(0 <= (cols <- as.integer(cols)), cols < n)
42      stopifnot(is.character(shape), nchar(shape) == 1,      m <- length(cols)
43                any(shape == c("t","s","g"))) # triangular / symmetric / general      if(missing(kind))
44      kind <-      kind <-
45          if(is.double(x)) "d"          if(is.double(x)) "d"
46          else if(is.logical(x)) "l"          else if(is.logical(x)) "l"
# Line 45  Line 48 
48              storage.mode(x) <- "double"              storage.mode(x) <- "double"
49              "d"              "d"
50          }          }
51      new(paste(kind, shape, "CMatrix", sep=''),      else stopifnot(any(kind == c("d","l","n")))
52          Dim = c(n,n), x = x, uplo = uplo,      if(kind != "n") {
53          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)          if((lx <- length(x)) == 1) x <- rep.int(x, m)
54            else if(lx != m) stop("length(x) must be either 1 or #{cols}")
55        }
56        stopifnot(is.character(shape), nchar(shape) == 1,
57                  any(shape == c("t","s","g"))) # triangular / symmetric / general
58        if(kind == "n") {
59            if(shape == "g")
60                new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
61            else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
62                     i = cols, p = 0:m)
63        }
64        ## kind != "n" -- have x slot :
65        else if(shape == "g")
66            new(paste0(kind, "gCMatrix"), Dim = c(n,m),
67                x = x, i = cols, p = 0:m)
68        else new(paste0(kind, shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
69                 x = x, i = cols, p = 0:m)
70  }  }
71    
72  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
73  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
74      .sparseDiagonal(n, x, uplo, shape = "s")      .sparseDiagonal(n, x, uplo, shape = "s")
75    
76  ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:  # instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
77  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")  .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
78      .sparseDiagonal(n, x, uplo, shape = "t")      .sparseDiagonal(n, x, uplo, shape = "t")
79    
# Line 217  Line 236 
236  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ddiMatrix", "sparseMatrix", diag2tT)
237  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
238  setAs("ddiMatrix", "TsparseMatrix", diag2tT)  setAs("ddiMatrix", "TsparseMatrix", diag2tT)
239    setAs("ddiMatrix", "dsparseMatrix", diag2tT)
240  setAs("ddiMatrix", "CsparseMatrix",  setAs("ddiMatrix", "CsparseMatrix",
241        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "d"), "CsparseMatrix"))
242  setAs("ddiMatrix", "symmetricMatrix",  setAs("ddiMatrix", "symmetricMatrix",
# Line 228  Line 248 
248  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)  ##_no_longer_ setAs("ldiMatrix", "sparseMatrix", diag2tT)
249  ## needed too (otherwise <dense> -> Tsparse is taken):  ## needed too (otherwise <dense> -> Tsparse is taken):
250  setAs("ldiMatrix", "TsparseMatrix", diag2tT)  setAs("ldiMatrix", "TsparseMatrix", diag2tT)
251    setAs("ldiMatrix", "lsparseMatrix", diag2tT)
252  setAs("ldiMatrix", "CsparseMatrix",  setAs("ldiMatrix", "CsparseMatrix",
253        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))        function(from) as(.diag2tT(from, "U", "l"), "CsparseMatrix"))
254  setAs("ldiMatrix", "symmetricMatrix",  setAs("ldiMatrix", "symmetricMatrix",
# Line 334  Line 355 
355            function(x = 1, nrow, ncol) .diag.x(x))            function(x = 1, nrow, ncol) .diag.x(x))
356    
357  subDiag <- function(x, i, j, ..., drop) {  subDiag <- function(x, i, j, ..., drop) {
358      x <- as(x, "TsparseMatrix")      x <- as(x, "CsparseMatrix") ## << was "TsparseMatrix" (Csparse is faster now)
359      x <- if(missing(i))      x <- if(missing(i))
360          x[, j, drop=drop]          x[, j, drop=drop]
361      else if(missing(j))      else if(missing(j))
# Line 602  Line 623 
623      ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])      ind <- rep.int(seq_len(dx[2]), x@p[-1] - x@p[-dx[2]-1L])
624      if(y@diag == "N")      if(y@diag == "N")
625          x@x <- x@x * y@x[ind]          x@x <- x@x * y@x[ind]
626        if(is(x, "compMatrix") && length(xf <- x@factors)) {
627            ## instead of dropping all factors, be smart about some
628            ## TODO ......
629            x@factors <- list()
630        }
631      x      x
632  }  }
633    
# Line 611  Line 637 
637      if(dx[2] != dy[1]) stop("non-matching dimensions")      if(dx[2] != dy[1]) stop("non-matching dimensions")
638      if(x@diag == "N")      if(x@diag == "N")
639          y@x <- y@x * x@x[y@i + 1L]          y@x <- y@x * x@x[y@i + 1L]
640        if(is(y, "compMatrix") && length(yf <- y@factors)) {
641            ## instead of dropping all factors, be smart about some
642            ## TODO
643            keep <- character()
644            if(iLU <- names(yf) == "LU") {
645                ## TODO keep <- "LU"
646            }
647            y@factors <- yf[keep]
648        }
649      y      y
650  }  }
651    
# Line 704  Line 739 
739          .diag.2N(e1)          .diag.2N(e1)
740      }      }
741      else { ## result not diagonal, but at least symmetric:      else { ## result not diagonal, but at least symmetric:
742            ## e.g., m == m
743          isNum <- (is.numeric(r) || is.numeric(r00))          isNum <- (is.numeric(r) || is.numeric(r00))
744          isLog <- (is.logical(r) || is.logical(r00))          isLog <- (is.logical(r) || is.logical(r00))
745            Matrix.msg("exploding <diag> o <diag> into dense matrix", .M.level = 2)
         Matrix.msg("exploding   <diag>  o  <diag>  into dense matrix")  
746          d <- e1@Dim          d <- e1@Dim
747          n <- d[1]          n <- d[1]
748          stopifnot(length(r) == n)          stopifnot(length(r) == n)
# Line 735  Line 770 
770              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)              setMethod("Ops", signature(e1 = c1, e2 = c2), diagOdiag)
771  }  }
772    
773  ## FIXME:    diagonal  o  triangular  |-->  triangular  ## diagonal  o  triangular  |-->  triangular
774  ## -----     diagonal  o  symmetric   |-->  symmetric  ## diagonal  o  symmetric   |-->  symmetric
775  ##    {also when other is sparse: do these "here" --  ##    {also when other is sparse: do these "here" --
776  ##     before conversion to sparse, since that loses "diagonality"}  ##     before conversion to sparse, since that loses "diagonality"}
777    
# Line 921  Line 956 
956  for(other in c("ANY", "Matrix", "dMatrix")) {  for(other in c("ANY", "Matrix", "dMatrix")) {
957      ## ddi*:      ## ddi*:
958      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ddiMatrix", e2 = other),
959                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "d"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "d"), e2))
960      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ddiMatrix"),
961                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "d")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "d")))
962      ## ldi*:      ## ldi*:
963      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),      setMethod("Ops", signature(e1 = "ldiMatrix", e2 = other),
964                function(e1,e2) callGeneric(diag2tT.u(e1,e2, "l"), e2))                function(e1,e2) callGeneric(diag2Tsmart(e1,e2, "l"), e2))
965      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),      setMethod("Ops", signature(e1 = other, e2 = "ldiMatrix"),
966                function(e1,e2) callGeneric(e1, diag2tT.u(e2,e1, "l")))                function(e1,e2) callGeneric(e1, diag2Tsmart(e2,e1, "l")))
967  }  }
968    
969  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :  ## Direct subclasses of "denseMatrix": currently ddenseMatrix, ldense... :
970    if(FALSE)## too general, would contain  denseModelMatrix:
971  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses  dense.subCl <- local({ dM.scl <- getClass("denseMatrix")@subclasses
972                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })                         names(dM.scl)[sapply(dM.scl, slot, "distance") == 1] })
973    dense.subCl <- paste(c("d","l","n"), "denseMatrix", sep="")
974  for(DI in diCls) {  for(DI in diCls) {
975      for(c2 in c(dense.subCl, "Matrix")) {      for(c2 in c(dense.subCl, "Matrix")) {
976          for(Fun in c("*", "^", "&")) {          for(Fun in c("*", "^", "&")) {

Legend:
Removed from v.2458  
changed lines
  Added in v.2608

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