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 2482, Wed Sep 30 09:37:15 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      ii <- if(n) 0:(n - 1L) else integer(0)      else stopifnot(any(kind == c("d","l","n")))
52        if(kind != "n") {
53            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")      if(shape == "g")
60          new(paste0(kind, "gCMatrix"), Dim = c(n,n),              new("ngCMatrix", Dim = c(n,m), i = cols, p = 0:m)
61              x = x, i = ii, p = 0:n)          else new(paste0("n", shape, "CMatrix"), Dim = c(n,m), uplo = uplo,
62      else new(paste0(kind, shape, "CMatrix"), Dim = c(n,n), uplo = uplo,                   i = cols, p = 0:m)
63               x = x, i = ii, p = 0:n)      }
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 339  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 607  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 616  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 709  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 740  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 926  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.2482  
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