spdiags grinds to a halt!

I'm trying to insert 441 diagonals into a sparse matrix that is 614292x614292 in size. My matrix vaguely a Toeplitz matrix and I have something like 32 GB of ram.

After inserting about 150 diagonals the inserting slows to a halt.

Looking at spdiags it appears to insert each diagonal seperately and perform some kind of 'compression', I'm wondering if there is some way to insert them all at once.

function S=buildconvmatrix(Cb,m,n)
%In the future have custom diagonals...
[cm,cn]=size(Cb);
[Xi,Yi]=meshgrid(1:cm,1:cn);
Xi=reshape(Xi,[],1);
Yi=reshape(Yi,[],1);
zpoint = ceil(numel(Cb)/2);
S = sparse(m*n,m*n);
for q=1:numel(Xi)
    row(q) = q-zpoint;
end
vList = reshape(Cb,1,[]);
S=spdiags(ones(m*n,1)*vList,row,S);% Hours
end

--Edit--

For example this is much faster, which leads to me believe that there is some kind of performance bug in Matlab's stock implementation.

function S=buildconvmatrixfaster(Cb,m,n)
[cm,cn]=size(Cb);
[Xi,Yi]=meshgrid(1:cm,1:cn);
Xi=reshape(Xi,[],1);
Yi=reshape(Yi,[],1);
zpoint = ceil(numel(Cb)/2);
cuts = 12;
M=cell(cuts,1);
for i=1:cuts
    M{i}=sparse(m*n,m*n);
end
els=numel(Xi);
bin=els/cuts;
for q=1:els
    disp(strcat('On:',num2str(q),'/',num2str(numel(Xi))));
    v=Cb(Xi(q),Yi(q));
    row = q-zpoint;
    foo=ceil(q/bin);
    M{foo}=spdiags(ones(m*n,1)*v,row,M{foo});
end
S=sparse(m*n,m*n);
disp('Adding');
for i=1:cuts
S=S+M{i};
end
end

[Mikhail Kandel]

Also there is probably an opportunity for parallelism here because you can make two sparse matrices and add them?

[Matt J]

What are the dimensions of Cb?

[Mikhail Kandel]

@Matt J its 21x21

Answer 1

[Matt J]

Instead of

S = sparse(m*n,m*n);

try

S = spalloc(m*n,m*n, numel(ones(m*n,1)*vList));

[Mikhail Kandel]

This had no measurable effect.

Answer 2

[Matt J]

See also if this

http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions

is applicable to you.

[Mikhail Kandel]

This might be close, I'm trying to check a C++ program where I hit each pixel of an image with a varying convolution kernel. I'm writing it in matlab to get the more natural mathematical form, so I can debug if my math is wrong or if my C++ is wrong.