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 2207, Mon Jul 7 22:34:52 2008 UTC revision 2237, Fri Jul 25 06:55:42 2008 UTC
# Line 32  Line 32 
32      }      }
33  }  }
34    
35  ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()  .sparseDiagonal <- function(n, x = rep.int(1,n), uplo = "U", shape = "t") {
 .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U") {  
36      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)      stopifnot(n == (n. <- as.integer(n)), (n <- n.) >= 0)
37      if((lx <- length(x)) == 1) x <- rep.int(x, n)      if((lx <- length(x)) == 1) x <- rep.int(x, n)
38      else if(lx != n) stop("length(x) must be 1 or n")      else if(lx != n) stop("length(x) must be 1 or n")
39      cls <-      stopifnot(is.character(shape), nchar(shape) == 1,
40          if(is.double(x)) "dsCMatrix"                any(shape == c("t","s","g"))) # triangular / symmetric / general
41          else if(is.logical(x)) "lsCMatrix"      kind <-
42            if(is.double(x)) "d"
43            else if(is.logical(x)) "l"
44          else { ## for now          else { ## for now
45              storage.mode(x) <- "double"              storage.mode(x) <- "double"
46              "dsCMatrix"              "d"
47          }          }
48      new(cls, Dim = c(n,n), x = x, uplo = uplo,      new(paste(kind, shape, "CMatrix", sep=''),
49            Dim = c(n,n), x = x, uplo = uplo,
50          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)          i = if(n) 0:(n - 1L) else integer(0), p = 0:n)
51  }  }
52    
53    ## Pkg 'spdep' had (relatively slow) versions of this as_dsCMatrix_I()
54    .symDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
55        .sparseDiagonal(n, x, uplo, shape = "s")
56    
57    ## instead of   diagU2N(as(Diagonal(n), "CsparseMatrix")), diag = "N" in any case:
58    .trDiagonal <- function(n, x = rep.int(1,n), uplo = "U")
59        .sparseDiagonal(n, x, uplo, shape = "t")
60    
61    
62  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.  ### This is modified from a post of Bert Gunter to R-help on  1 Sep 2005.
63  ### Bert's code built on a post by Andy Liaw who most probably was influenced  ### Bert's code built on a post by Andy Liaw who most probably was influenced
64  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002  ### by earlier posts, notably one by Scott Chasalow on S-news, 16 Jan 2002
65  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
66  if(FALSE)##--- currently unused:  if(FALSE)##--- no longer used:
67  .bdiag <- function(lst) {  .bdiag <- function(lst) {
68      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ### block-diagonal matrix [a dgTMatrix] from list of matrices
69      stopifnot(is.list(lst), length(lst) >= 1)      stopifnot(is.list(lst), length(lst) >= 1)
# Line 74  Line 85 
85      }      }
86      r      r
87  }  }
88  ### Doug Bates needed something like bdiag() for lower-triangular  ## expand(<mer>) needed something like bdiag() for lower-triangular
89  ### (Tsparse) Matrices and provided a much more efficient implementation:  ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
90    ##  implementation for those; now extended and generalized:
91  .bdiag <- function(lst) {  .bdiag <- function(lst) {
92      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
93      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
94    
95      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
96                     as, "TsparseMatrix")                     as, "TsparseMatrix")
   
97      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
98      ## else      ## else
99      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
100      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
101      new("dgTMatrix", Dim = c(i_off[nl+1], j_off[nl + 1]),  
102          i = unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k])),      clss <- sapply(Tlst, class)
103          j = unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])),      knds <- substr(clss, 2, 2)
104          x = unlist(lapply(Tlst, slot, "x")))      sym  <- knds == "s" # symmetric ones
105        tri  <- knds == "t" # triangular ones
106        use.n <- any(is.n <- substr(clss,1,1) == "n")
107        if(use.n && !(use.n <- all(is.n)))
108            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
109        if(all(sym)) { ## result should be *symmetric*
110            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
111            tLU <- table(uplos)# of length 1 or 2 ..
112            if(length(tLU) == 1) { ## all "U" or all "L"
113                useU <- uplos[1] == "U"
114            } else { ## length(tLU) == 2, counting "L" and "U"
115                useU <- diff(tLU) >= 0
116                if(useU && (hasL <- tLU[1] > 0))
117                    Tlst[hasL] <- lapply(Tlst[hasL], t)
118                else if(!useU && (hasU <- tLU[2] > 0))
119                    Tlst[hasU] <- lapply(Tlst[hasU], t)
120            }
121            if(use.n) { ## return nsparseMatrix :
122                r <- new("nsTMatrix")
123            } else {
124                r <- new("dsTMatrix")
125                r@x <- unlist(lapply(Tlst, slot, "x"))
126            }
127            r@uplo <- if(useU) "U" else "L"
128        }
129        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
130                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
131           ){ ## *triangular* result
132    
133            if(use.n) { ## return nsparseMatrix :
134                r <- new("ntTMatrix")
135            } else {
136                r <- new("dtTMatrix")
137                r@x <- unlist(lapply(Tlst, slot, "x"))
138            }
139            r@uplo <- ULs[1L]
140        }
141        else {
142            if(any(sym))
143                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
144            if(use.n) { ## return nsparseMatrix :
145                r <- new("ngTMatrix")
146            } else {
147                r <- new("dgTMatrix")
148                r@x <- unlist(lapply(Tlst, slot, "x"))
149            }
150        }
151        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
152        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
153        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
154        r
155  }  }
156    
157  bdiag <- function(...) {  bdiag <- function(...) {
# Line 181  Line 242 
242                Dim = from@Dim, Dimnames = from@Dimnames)                Dim = from@Dim, Dimnames = from@Dimnames)
243        })        })
244    
245    setAs("diagonalMatrix", "nsparseMatrix", function(from) as(from, "nMatrix"))
246    
247  ## Cheap fast substitute for diag() which *does* preserve the mode of x :  ## Cheap fast substitute for diag() which *does* preserve the mode of x :
248  mkDiag <- function(x, n) {  mkDiag <- function(x, n) {
# Line 899  Line 961 
961      invisible(x)      invisible(x)
962  }  }
963    
964    ## somewhat consistent with "print" for sparseMatrix :
965    setMethod("print", signature(x = "diagonalMatrix"), prDiag)
966    
967  setMethod("show", signature(object = "diagonalMatrix"),  setMethod("show", signature(object = "diagonalMatrix"),
968            function(object) {            function(object) {
969                d <- dim(object)                d <- dim(object)

Legend:
Removed from v.2207  
changed lines
  Added in v.2237

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