SCM Repository

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

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

revision 2203, Sat Jun 14 20:09:17 2008 UTC revision 2206, Wed Jun 25 14:57:02 2008 UTC
# Line 52  Line 52
52  ### 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
53  ### 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
54  ### who posted his bdiag() function written in December 1995.  ### who posted his bdiag() function written in December 1995.
55    if(FALSE)##--- currently unused:
56  bdiag <- function(...) {  .bdiag <- function(lst) {
57      if(nargs() == 0) return(new("dgCMatrix"))      ### block-diagonal matrix [a dgTMatrix] from list of matrices
58      ## else :      stopifnot(is.list(lst), length(lst) >= 1)
59      mlist <- if (nargs() == 1) as.list(...) else list(...)      dims <- sapply(lst, dim, USE.NAMES=FALSE)
dims <- sapply(mlist, dim)
60      ## make sure we had all matrices:      ## make sure we had all matrices:
61      if(!(is.matrix(dims) && nrow(dims) == 2))      if(!(is.matrix(dims) && nrow(dims) == 2))
62          stop("some arguments are not matrices")          stop("some arguments are not matrices")
63      csdim <- rbind(rep.int(0L, 2),      csdim <- rbind(rep.int(0L, 2),
64                     apply(sapply(mlist, dim), 1, cumsum))                     apply(dims, 1, cumsum))
65      ret <- new("dgTMatrix", Dim = as.integer(csdim[nrow(csdim),]))      r <- new("dgTMatrix")
66        r@Dim <- as.integer(csdim[nrow(csdim),])
67      add1 <- matrix(1:0, 2,2)      add1 <- matrix(1:0, 2,2)
68      for(i in seq_along(mlist)) {      for(i in seq_along(lst)) {
69          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])          indx <- apply(csdim[i:(i+1),] + add1, 2, function(n) n[1]:n[2])
70          if(is.null(dim(indx))) ## non-square matrix          if(is.null(dim(indx))) ## non-square matrix
71              ret[indx[[1]],indx[[2]]] <- mlist[[i]]              r[indx[[1]],indx[[2]]] <- lst[[i]]
72          else ## square matrix          else ## square matrix
73              ret[indx[,1],indx[,2]] <- mlist[[i]]              r[indx[,1], indx[,2]] <- lst[[i]]
74        }
75        r
76      }      }
77      ## slightly debatable if we really should return Csparse.. :  ### Doug Bates needed something like bdiag() for lower-triangular
78      as(ret, "CsparseMatrix")  ### (Tsparse) Matrices and provided a much more efficient implementation:
79    .bdiag <- function(lst) {
80        ### block-diagonal matrix [a dgTMatrix] from list of matrices
81        stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
82
83        Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
84                       as, "TsparseMatrix")
85
86        if(nl == 1) return(Tlst[[1]])
87        ## else
88        i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
89        j_off <- c(0L, cumsum(sapply(Tlst, ncol)))
90        new("dgTMatrix", Dim = c(i_off[nl+1], j_off[nl + 1]),
91            i = unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k])),
92            j = unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])),
93            x = unlist(lapply(Tlst, slot, "x")))
94    }
95
96    bdiag <- function(...) {
97        if((nA <- nargs()) == 0) return(new("dgCMatrix"))
98        if(nA == 1 && !is.list(...))
99            return(as(..., "CsparseMatrix"))
100        alis <- if(nA == 1 && is.list(..1)) ..1 else list(...)
101        if(length(alis) == 1)
102            return(as(alis[[1]], "CsparseMatrix"))
103
104        ## else : two or more arguments
105        as(.bdiag(alis), "CsparseMatrix")
106  }  }
107
108

Legend:
 Removed from v.2203 changed lines Added in v.2206

 root@r-forge.r-project.org ViewVC Help Powered by ViewVC 1.0.0
Thanks to: