dan
dan

Reputation: 41

Optimal way of doing iterative assembly of sparse matrices in Matlab?

My code needs to in a loop modify the elements of a sparse matrix. Doing this matlab warns me that This sparse indexing expression is likely to be slow. I am preallocating the sparse arrays using the Spalloc function but am still getting this warning. What is the optimal approach for assembling of sparse matrices? This is what I am currently doing.

K=spalloc(n,n,100); f=spalloc(n,1,100);

for i = 1:Nel
  [Ke,fe] = myFunction(Ex(i),Ey(i));
  inds = data(i,2:end);
  K(inds,inds) = K(inds,inds) + Ke; 
  f(inds) = f(inds)+fe; 
end

the indices in inds may be appear several times in the loop, meaning an element in K or f may receive multiple contributions. The last two lines within the loop are where I'm getting warnings.

Upvotes: 0

Views: 241

Answers (1)

Ryan Livingston
Ryan Livingston

Reputation: 1928

A common approach is to use the triplet form of the sparse constructor:

S = sparse(i,j,v,m,n)

i and j are row and column index vectors and v is the corresponding data vector. Values corresponding to repeated indices are summed like your code does. So you could instead build up row and column index vectors along with a data vector and then just call sparse with those.

For example something like:

nout = Nel*(size(data,2)-1);
% Data vector for K
Kdata = zeros(1,nout);
% Data vector for f
fdata = zeros(1,nout);
% Index vector for K and f
sparseIdxvec = ones(1,nout);
outIdx = 1;
for i = 1:Nel
  [Ke,fe] = myFunction(Ex(i),Ey(i));
  inds = data(i,2:end);
  nidx = numel(inds);
  outIdxvec = outIdx:outIdx+nidx-1;
  sparseIdxvec(outIdxvec) = inds;
  Kdata(outIdxvec) = Ke;
  fdata(outIdxvec) = fe;
  outIdx = outIdx + nidx;
end
K = sparse(sparseIdxvec,sparseIdxvec,Kdata,n,n);
f = sparse(sparseIdxvec,1,fdata,n,1);

Depending on your application, that may or may not actually be faster.

Upvotes: 0

Related Questions