matlab: extract block diagonals of large sparse ma

I have a large sparse matrix A, and I would like to create a sparse matrix of the 3X3 block diagonals of A. How would I do this? keep in mind that A is very large and sparse, so any methods that use iteration will be slow, and any methods that use some methods that creates full (as opposed to sparse) matrices will take up too much memory.

If I understand correctly, here is some code (see the portions between the %%%%%%%%%%% lines. Below are timing results, which seem reasonable to me, despite the for loop. The only trick is the use of the spalloc function, which you may have to tune for your application.

for N= [(3:3:12) (15:600:9000)]    
    bigsparse = sprand(N,N,0.1);
    tic;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    origSize = size(bigsparse);
    diagSize = 3;
    numDiags = size(bigsparse,1)/diagSize;
    assert(numDiags == floor(numDiags))

    bigsparse_diagonals = spalloc(origSize(1), origSize(2), ceil(prod(origSize)*0.1));
    for ix=(1:numDiags)-1
        ixsCurrent = ix*diagSize+[1:diagSize];
        bigsparse_diagonals(ixsCurrent,ixsCurrent) = ...
            bigsparse(ixsCurrent,ixsCurrent);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    fprintf(1,'%5d size --> %6.5f seconds \n', N, toc)
end

Timing results (note, it actually takes a lot longer to generate the random test matrix than to do the reformatting):

    3 size --> 0.00135 seconds 
    6 size --> 0.00014 seconds 
    9 size --> 0.00013 seconds 
   12 size --> 0.00014 seconds 
   15 size --> 0.00015 seconds 
  615 size --> 0.00392 seconds 
 1215 size --> 0.00874 seconds 
 1815 size --> 0.01537 seconds 
 2415 size --> 0.02570 seconds 
 3015 size --> 0.03595 seconds 
 3615 size --> 0.05007 seconds 
 4215 size --> 0.06420 seconds 
 4815 size --> 0.08690 seconds 
 5415 size --> 0.10077 seconds 
 6015 size --> 0.13322 seconds 
 6615 size --> 0.14923 seconds 
 7215 size --> 0.17562 seconds 
 7815 size --> 0.37371 seconds 
 8415 size --> 0.23060 seconds