Folks: In answer to a query, Andy Liaw recently submitted some code to construct a block diagonal matrix. For what seemed a fairly straightforward task, the code seemed a little "overweight" to me (that's an American stock analyst's term, btw), so I came up with a slightly cleaner version (with help from Andy):
bdiag<-function(...){ mlist<-list(...) ## handle case in which list of matrices is given if(length(mlist)==1)mlist<-unlist(mlist,rec=FALSE) csdim<-rbind(c(0,0),apply(sapply(mlist,dim),1,cumsum )) ret<-array(0,dim=csdim[length(mlist)+1,]) add1<-matrix(rep(1:0,2),nc=2) for(i in seq(along=mlist)){ indx<-apply(csdim[i:(i+1),]+add1,2,function(x)x[1]:x[2]) ## non-square matrix if(is.null(dim(indx)))ret[indx[[1]],indx[[2]]]<-mlist[[i]] ## square matrix else ret[indx[,1],indx[,2]]<-mlist[[i]] } ret } I doubt that there's any noticeable practical performance difference, of course. The strategy is entirely basic: just get the right indices for replacement of the arguments into a matrix of 0's of the right dimensions. About the only thing to notice is that the apply() construction returns either a list or matrix depending on whether a matrix is square or not (a subtlety that tripped me up in my first version of this code). I would be pleased if this stimulated others to come up with cleverer/more elegant approaches that they would share, as it's the sort of thing that I'll learn from and find useful. Cheers to all, Bert Gunter ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html