Sparse FEM Matrix build inefficiency

Hello everyone.

I am trying to improve the efficiency of my FEM matlab code using sparse matrices. I build three vectors sp_i sp_j and sp_data and then call the sparse(sp_i,sp_j,sp_data,dofs,dofs) ONLY once to build the sparse matrix. However, the problem that occurs is that the build of the sp_data vector is inefficient and the developed code is much slower than the original one (with the dense matrix inversion).

Here is the code (two versions):

(practically A: the assembly is performed adding a 3by3 matrix to the dense dofs_by_dofs global matrix each time B: the assembly is performed adding 3 3by1 vectors at the global sp_data vector each time )

A) Original assembly (dense):

for h1=1:1:(hexa1)
for h2=1:1:(hexa2)
for h3=1:1:(hexa3)
elem=h1+(h2-1)*hexa1+(h3-1)*(hexa1)*hexa2; 
for line=1:1:8
for row=1:1:8
K_sunol_disp(((t_disp(elem,line)-1)*3+1):((t_disp(elem,line)-1)*3+3),((t_disp(elem,row)-1)*3+1):((t_disp(elem,row)-1)*3+3))=...
K_sunol_disp(((t_disp(elem,line)-1)*3+1):((t_disp(elem,line)-1)*3+3),((t_disp(elem,row)-1)*3+1):((t_disp(elem,row)-1)*3+3))+...
k_stoixeiou(((line-1)*3+1):((line-1)*3+3),((row-1)*3+1):((row-1)*3+3),elem);
end        
end

end
end
end
t_disp is a matrix that contains connectivity data of the FEM degrees of freedom.

B)Code development (assembly of sparse(...) vectors):

for h1=1:1:(hexa1)
for h2=1:1:(hexa2)
for h3=1:1:(hexa3)
elem=h1+(h2-1)*hexa1+(h3-1)*(hexa1)*hexa2;
for line=1:1:8
for row=1:1:8
thesi_dedom=t_K_sunol(komvoi_8+t_disp(elem,line),komvoi_8+t_disp(elem,row));
line_31=(line-1)*3+1;
row_3=(row-1)*3;
sp_data(thesi_dedom:thesi_dedom+8,1)=sp_data(thesi_dedom:thesi_dedom+8,1)+[k_stoixeiou((line_31):(line_31+2),(row_3+1),elem);...
                                k_stoixeiou((line_31):(line_31+2),(row_3+2),elem);k_stoixeiou((line_31):(line_31+2),(row_3+3),elem)];
end        
end

end
end
end
K_sunol_disp=sparse(sp_i,sp_j,sp_data,dofs,dofs);

The most costly command as indicated by the "run and time" is the "sp_data(thesi_dedom:thesi_dedom+8,1)=sp_data(thesi_dedom:thesi_dedom+8,1)+....". Dividing it in three simple commands didn't help at all. sp_data is of course preallocated.

Is there a command that has as input the position vectors sp_i and sp_j and a stored as dense matrix K_sunol_disp (not fully populated) that will be turned into sparse?

Is there a way to prevent the build of the sp_data vector from slowing down the code that will allow us to benefit from the sparsity of the matrices? (I will also like to say that sp_i and sp_j don't create any problem as they are built once and then used multiple times(every time that I assemble the sparse matrix in the iterative solution procedure)).

Is there a better way to structure the problem solution ?

I would really appreciate any advise because the addressing of this inefficiency would make the use of the FEM code that I develop feasible.

Thank you in advance

Answer 1

[Precise Simulation]

Although looping over all element, finite element forms, and quadrature points is the typical fem assembly approach, even with the Matlab JIT compilation this will be extremely inefficient. If you vectorize your assembly loop you can even approach the efficiency of Fortran (at the cost of using more memory) as done in this fem solver and assembly benchmark.