Optimized filling of sparse matrix using spalloc

Hi Everyone, A similar question has been asked on this forum but I think we need to revisit the answers provided. My objective is as follows:

%Create allocated memory for row, column, and value vectors
row_vec=spalloc(n*n,1,ceil(n*n*0.004));col_vec=spalloc(n*n,1,ceil(n*n*0.004)); 
val_vec=spalloc(n*n,1,ceil(n*n*0.004));
%Looping over the number of simplex elements, in this case linear tetrahedra
for ii=1:nt
row_vec(k,1)=p;col_vec(k,1)=q;val_vec(k,1)=val_vec(k,1)+r;
end
[rrnz,rcnz,rvnz]=find(row_vec);[crnz,ccnz,cvnz]=find(col_vec);[vrnz,vcnz,vvnz]=find(val_vec);
rowfinal=sparse(rrnz,rcnz,rvnz,nnd*nnd,1);colfinal=sparse(crnz,ccnz,cvnz,nnd*nnd,1);
vfinal=sparse(vrnz,vcnz,vvnz,nnd*nnd,1);
A_matrix=sparse(rowfinal,colfinal,vfinal,n,n);

******************************************************************************************

I see that row_vec,col_vec, and val_vec are still sparse matrices being filled in a "slow" manner. I have to derive the p,q, and r values within the loop. Is there a more efficient way to fill row_vec, col_vec, and val_vec within the loop? Many thanks. Souvik

Answer 1

[John D'Errico]

DON'T FILL A SPARSE MATRIX SLOWLY. You will virtually never be happy with the result. DON'T create those row and column arrays as sparse. Again, not a good idea. As you are collecting your matrices, this is the wrong way. A bad idea.

You have a known set of tetrahedra (nt of them.) So why do you need to accumulate the data using sparse matrices? Sorry, but that is just plain silly.

Even if you don't know the final size of your arrays, there are better methods. But here you know how big they will be.

It is difficult to know exactly what you are doing, because your code is incomplete. But this should give you the basic idea.

spdata = zeros(nt,3);
for ii=1:nt
  % stuff where you compute p,q,r
  spdata(ii,:) = [p,q,r];
end
A_matrix=sparse(spdata(:,1),spdata(:,2),spdata(:,3),n,n);

Note that sparse is smart enough to accumulate values that happen to lie at the same index. Thus if p and q are identical for two cases, then the corresponding values will be summed up.

As I said, in the event where nt is not in fact known in advance, there are ways to grow that array efficiently. But here you know it.

[Souvik Mukherjee]

John, You are really smart. However, the matrix size I need to fill is nXn where n is the number of all nodes in the tetrahedra. If I define 3 single column zero vectors that are (nxn,1) size, matlab runs out of memory. If I define zeros(nXn,3) I have the same issue. The best way out is to define the index and value vectors as sparse (nXn)x1 column vectors with max nonzero elements not exceeding 0.4% of the total number of elements. The loop over nt which is the the number of tetrahedra allows me to find out the value r for row p and col q for tetrahdera i. The total value for the (p,q) element of the A_matrix that I am assembling (outside the loop) is given by summing over all the tetrahedra where both p and q nodes are present. Hope this helps you understand the problem better. Thanks.