[R] More block diagonal matrix construction code

[Berton Gunter]

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