require(foreach)
require(itertools)

# block matrix multiplication
bmm <- function(x, y, chunkSize=max(1, floor(nrow(x)/6)), verbose=FALSE,
  writeBlock=FALSE, returnResult=TRUE, filePrefix="block.", writeDir="",
  projectionBlocks=FALSE) {
  if (ncol(x) != nrow(y))
    stop("Non-conformable matrices")
  if (verbose) {
    
    cat(ceiling(nrow(x)/chunkSize)^2, "block computations to perform.\n")
  }
  foreach(i=isplitVector(1:nrow(x), chunkSize=chunkSize), ic=icount(),
    .combine=rbind) %do% {
    ret <- foreach (j=isplitVector(1:ncol(y), chunkSize=chunkSize), 
      jc=icount(), .combine=cbind ) %dopar% {
      if (projectionBlocks) {
        blockMult <- foreach(k=isplitVector(1:ncol(x), chunkSize=chunkSize),
          .combine="+") %do% {
          x[i,k, drop=FALSE] %*% y[k,j, drop=FALSE]
        }
      } else {
        blockMult <- x[i,, drop=FALSE] %*% y[,j,drop=FALSE]
      }
      if (writeBlock) {
        saveRDS(blockMult, 
          paste(writeDir, filePrefix, ic, ".", jc, ".rds", sep=""))
      }
      if (!returnResult)
        blockMult <- 1
    }
    ret
  }
}