Function bdiag to build the canonical form of a matrix

bdiag.R

######################################################################
# Function bdiag to build the canonical form of a matrix
#---------------------------------------------------------------------
# Pedro Jordano. 18 Dic 2007.
#---------------------------------------------------------------------
bdiag <- function(x){ 
     if(!is.list(x)) stop("x not a list") 
     n <- length(x) 
     if(n==0) return(NULL) 
     x <- lapply(x, function(y) if(length(y)) as.matrix(y) else 
stop("Zero-length component in x")) 
     d <- array(unlist(lapply(x, dim)), c(2, n)) 
     rr <- d[1,] 
     cc <- d[2,] 
     rsum <- sum(rr) 
     csum <- sum(cc) 
     out <- array(0, c(rsum, csum)) 
     ind <- array(0, c(4, n)) 
     rcum <- cumsum(rr) 
     ccum <- cumsum(cc) 
     ind[1,-1] <- rcum[-n] 
     ind[2,] <- rcum 
     ind[3,-1] <- ccum[-n] 
     ind[4,] <- ccum 
     imat <- array(1:(rsum * csum), c(rsum, csum)) 
     iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2], 
(y[3]+1):y[4]], imat=imat) 
     iuse <- as.vector(unlist(iuse)) 
     out[iuse] <- unlist(x) 
     return(out) 
} 
#### 
mats <- list(matrix(1:20, 5, 4), matrix(1:12, 4, 3), matrix(1:25, 5, 
5)) 
bdiag(mats) 
bdiag <- function(x){ 
+      if(!is.list(x)) stop("x not a list") 
+      n <- length(x) 
+      if(n==0) return(NULL) 
+      x <- lapply(x, function(y) if(length(y)) as.matrix(y) else 
+ stop("Zero-length component in x")) 
+      d <- array(unlist(lapply(x, dim)), c(2, n)) 
+      rr <- d[1,] 
+      cc <- d[2,] 
+      rsum <- sum(rr) 
+      csum <- sum(cc) 
+      out <- array(0, c(rsum, csum)) 
+      ind <- array(0, c(4, n)) 
+      rcum <- cumsum(rr) 
+      ccum <- cumsum(cc) 
+      ind[1,-1] <- rcum[-n] 
+      ind[2,] <- rcum 
+      ind[3,-1] <- ccum[-n] 
+      ind[4,] <- ccum 
+      imat <- array(1:(rsum * csum), c(rsum, csum)) 
+      iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2], 
+ (y[3]+1):y[4]], imat=imat) 
+      iuse <- as.vector(unlist(iuse)) 
+      out[iuse] <- unlist(x) 
+      return(out) 
+ }
#---------------------------------------------------------------------
R> ######################### 
R> mats <- list(matrix(1:20, 5, 4), matrix(1:12, 4, 3), matrix(1:25, 5, 
+ 5)) 
R> bdiag(mats) 
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,]    1    6   11   16    0    0    0    0    0     0     0     0
 [2,]    2    7   12   17    0    0    0    0    0     0     0     0
 [3,]    3    8   13   18    0    0    0    0    0     0     0     0
 [4,]    4    9   14   19    0    0    0    0    0     0     0     0
 [5,]    5   10   15   20    0    0    0    0    0     0     0     0
 [6,]    0    0    0    0    1    5    9    0    0     0     0     0
 [7,]    0    0    0    0    2    6   10    0    0     0     0     0
 [8,]    0    0    0    0    3    7   11    0    0     0     0     0
 [9,]    0    0    0    0    4    8   12    0    0     0     0     0
[10,]    0    0    0    0    0    0    0    1    6    11    16    21
[11,]    0    0    0    0    0    0    0    2    7    12    17    22
[12,]    0    0    0    0    0    0    0    3    8    13    18    23
[13,]    0    0    0    0    0    0    0    4    9    14    19    24
[14,]    0    0    0    0    0    0    0    5   10    15    20    25
R> m<-matrix(1:20, 5, 4)
R> ml<-list(m,t(m))
R> bdiag(ml)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]    1    6   11   16    0    0    0    0    0
 [2,]    2    7   12   17    0    0    0    0    0
 [3,]    3    8   13   18    0    0    0    0    0
 [4,]    4    9   14   19    0    0    0    0    0
 [5,]    5   10   15   20    0    0    0    0    0
 [6,]    0    0    0    0    1    2    3    4    5
 [7,]    0    0    0    0    6    7    8    9   10
 [8,]    0    0    0    0   11   12   13   14   15
 [9,]    0    0    0    0   16   17   18   19   20