# SCM Repository

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

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

revision 2207, Mon Jul 7 22:34:52 2008 UTC revision 2209, Mon Jul 14 19:39:39 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:  if(FALSE)##--- no longer used:
56  .bdiag <- function(lst) {  .bdiag <- function(lst) {
57      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ### block-diagonal matrix [a dgTMatrix] from list of matrices
58      stopifnot(is.list(lst), length(lst) >= 1)      stopifnot(is.list(lst), length(lst) >= 1)
# Line 74  Line 74
74      }      }
75      r      r
76  }  }
77  ### Doug Bates needed something like bdiag() for lower-triangular  ## expand(<mer>) needed something like bdiag() for lower-triangular
78  ### (Tsparse) Matrices and provided a much more efficient implementation:  ## (Tsparse) Matrices; hence Doug Bates provided a much more efficient
79    ##  implementation for those; now extended and generalized:
80  .bdiag <- function(lst) {  .bdiag <- function(lst) {
81      ### block-diagonal matrix [a dgTMatrix] from list of matrices      ## block-diagonal matrix [a dgTMatrix] from list of matrices
82      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)      stopifnot(is.list(lst), (nl <- length(lst)) >= 1)
83
84      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"      Tlst <- lapply(lapply(lst, Matrix:::as_Csp2), # includes "diagU2N"
85                     as, "TsparseMatrix")                     as, "TsparseMatrix")

86      if(nl == 1) return(Tlst[[1]])      if(nl == 1) return(Tlst[[1]])
87      ## else      ## else
88      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))      i_off <- c(0L, cumsum(sapply(Tlst, nrow)))
89      j_off <- c(0L, cumsum(sapply(Tlst, ncol)))      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])),      clss <- sapply(Tlst, class)
92          j = unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k])),      knds <- substr(clss, 2, 2)
93          x = unlist(lapply(Tlst, slot, "x")))      sym  <- knds == "s" # symmetric ones
94        tri  <- knds == "t" # triangular ones
95        use.n <- any(is.n <- substr(clss,1,1) == "n")
96        if(use.n && !(use.n <- all(is.n)))
97            Tlst[is.n] <- lapply(Tlst[is.n], as, "lMatrix")
98        if(all(sym)) { ## result should be *symmetric*
99            uplos <- sapply(Tlst, slot, "uplo") ## either "U" or "L"
100            tLU <- table(uplos)# of length 1 or 2 ..
101            if(length(tLU) == 1) { ## all "U" or all "L"
102                useU <- uplos[1] == "U"
103            } else { ## length(tLU) == 2, counting "L" and "U"
104                useU <- diff(tLU) >= 0
105                if(useU && (hasL <- tLU[1] > 0))
106                    Tlst[hasL] <- lapply(Tlst[hasL], t)
107                else if(!useU && (hasU <- tLU[2] > 0))
108                    Tlst[hasU] <- lapply(Tlst[hasU], t)
109            }
110            if(use.n) { ## return nsparseMatrix :
111                r <- new("nsTMatrix")
112            } else {
113                r <- new("dsTMatrix")
114                r@x <- unlist(lapply(Tlst, slot, "x"))
115            }
116            r@uplo <- if(useU) "U" else "L"
117        }
118        else if(all(tri) && { ULs <- sapply(Tlst, slot, "uplo")##  "U" or "L"
119                              all(ULs[1L] == ULs[-1L]) } ## all upper or all lower
120           ){ ## *triangular* result
121
122            if(use.n) { ## return nsparseMatrix :
123                r <- new("ntTMatrix")
124            } else {
125                r <- new("dtTMatrix")
126                r@x <- unlist(lapply(Tlst, slot, "x"))
127            }
128            r@uplo <- ULs[1L]
129        }
130        else {
131            if(any(sym))
132                Tlst[sym] <- lapply(Tlst[sym], as, "generalMatrix")
133            if(use.n) { ## return nsparseMatrix :
134                r <- new("ngTMatrix")
135            } else {
136                r <- new("dgTMatrix")
137                r@x <- unlist(lapply(Tlst, slot, "x"))
138            }
139        }
140        r@Dim <- c(i_off[nl+1], j_off[nl + 1])
141        r@i <- unlist(lapply(1:nl, function(k) Tlst[[k]]@i + i_off[k]))
142        r@j <- unlist(lapply(1:nl, function(k) Tlst[[k]]@j + j_off[k]))
143        r
144  }  }
145
146  bdiag <- function(...) {  bdiag <- function(...) {

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