# Efficient assembly of finite element matrix in MATLAB

Question

What is the most efficient algorithm for finding a row of a matrix which matches a given row? This is the same as a table lookup based on multiple criteria.

Context

Finite Element Matrices are usually large and sparse, so storing the entire matrix in memory is an inefficient use of computer memory. As such, I am generating the matrices as 3-tuples i.e. (row index, column index, matrix value). Due to the unstructured nature of the underlying grid, the same matrix element may need to be incremented a number of times.

To achieve this, the existing elements must be searched to see if the (row,column) pair already exists. If it does, then the value can be incremented. If it doesn't, a new 3-tuple is created for the new non-zero value of the matrix.

My question is: what is an efficient method for checking if the (row,column) pair already exists? And if it does, what the index into the array is?

This is essentially the same as finding a row of a matrix which matches a given row.

Currently, this is implemented in MATLAB as

% rw is the array of row indices
% cl is the array of column indices
% rv is the row being searched for
% cv is the column being searched for
% i.e. searching for (rv,cv) in (rw,cl)

possible_row = find(rw == rv);
column = find(cl(possible_row) == cv);

% check that it was found first
if isempty(column)
% doesn't exist yet
element_index = -1; % flag that it needs to be created
else
element_index = possible_row(column);
end

% return element_index


Due to the nature of find, this procedure is much faster than manually looping through the array. I could write a MEX function to do this, but I would prefer not to at this stage.

Any help will be greatly appreciated.

• This might be more appropriate for stackoverflow. Feb 2 '13 at 12:55
• @Alex Thanks. I'll leave it a day or so too see if there is any interest first. Feb 2 '13 at 22:40
• define large.... do you want an algorithm in MPI or just on a single core? Feb 6 '13 at 15:57
• Daryl, are you aware of the sparse matrix capability in MATLAB? If you truly want efficiency you should be using MATLAB sparse matrices, not your own tuple-based variant. Feb 6 '13 at 16:26
• To add to Aron's comment, MATLAB's sparse(row,col,val) will already automatically sum values for identical (row,col) pairs. Feb 6 '13 at 16:50

If you are assembling finite element matrices from arrays of row indices, column indices and values using MATLAB's sparse, there is no need to implement your own function to take care of repeated indices - sparse will automatically sum values with identical row and column indices (see the documentation, especially the third paragraph). Internally, this command creates a sparse compressed column matrix with the requested sparsity pattern and the given entries. For those interested in more details, this article describes the original implementation of MATLAB's sparse matrices.

This is indeed the fastest way of creating a sparse matrix for which you don't know the band structure (you can use spdiags for those); creating an empty matrix and assigning entries is much slower, since MATLAB needs to recreate the entire data structure after each assignment. Here's a small example to illustrate this (spalloc(n,n,3*n) allocates an empty sparse matrix with storage for 3n non-zero elements):

n = 20000;
%% Allocate, fill row-wise
tic;
A = spalloc(n,n,3*n);
for i = 1:n
A(i,1:3) = rand(3,1);
end
toc


Elapsed time is 6.738242 seconds.

%% Allocate, fill column-wise
tic;
A = spalloc(n,n,3*n);
for i = 1:n
A(1:3,i) = rand(3,1);
end
toc


Elapsed time is 1.524576 seconds.

%% Using sparse
tic;
row = zeros(3*n,1); col = zeros(3*n,1); val = zeros(3*n,1);
ind = 1;
for i = 1:n
row(ind+(0:2)) = [1,2,3];
col(ind+(0:2)) = [i,i,i];
val(ind+(0:2)) = rand(3,1);
ind = ind+3;
end
B = sparse(row,col,val,n,n);
toc


Elapsed time is 0.170769 seconds.

If you repeat that with n=40000, you'll see that the first two methods scale quadratically while the last method scales linearly.

• Thanks. This was very helpful. I wasn't aware of that functionality of sparse. Feb 7 '13 at 12:47

it might be a bit late, but we have some Matlab packages for it including 'classical' node elements and even for edge elements, see

http://www.mathworks.com/matlabcentral/fileexchange/27826-fast-assembly-of-stiffness-and-matrices-in-finite-element-method-using-nodal-elements

http://www.mathworks.com/matlabcentral/fileexchange/46635-fast-assembly-of-stiffness-and-matrices-in-finite-element-method-using-edge-elements

There are also reports enclosed, so implementation details should be clearer. The main message is, that these FEM assemblies can be done very fast in Matlab as well. For instance, you can easily generate a FEM matrix for 1 milion x 1 milion elements in several seconds, depending on your CPU speed.

Best, Jan

You are currently using a very inefficient data structure and algorithm. A simple analysis of your algorithm shows that it runs in $O(n_{nz}^2)$ time in the worst and average case, and since the find function iterates from the beginning of your array, if you are doing the standard forward loop over your basis elements, you are in the worst case scenario.

The most efficient general strategy for what you are trying to do in MATLAB (assemble a matrix from finite element basis) is to preallocate an empty sparse matrix using spalloc using whatever knowledge of the matrix's structure you have, then loop over the basis column-wise to assemble the matrix and simply use the the native MATLAB matrix syntax to update.

If you were working in C, C++, or Fortran, I would have some other advice about how to make your code cache friendly, but that really isn't a detail most MATLAB developers should be paying attention to.

• I'm sorry, but I have to disagree here: Assigning previously zero elements is pretty much the most expensive way of building a sparse matrix in MATLAB, because that forces a rebuild of the sparsity structure of the matrix in every iteration. Much faster in general is filling arrays of row indices, column indices and values in a (single) loop and then calling sparse. (In fact, this scales linearly in the number of columns, while your way scales quadratically.) But your point of letting MATLAB worry about the low level stuff is still valid, of course. Feb 6 '13 at 22:52
• From the implementations that I have used prior to the current version, this method was not efficient at all. The sparsity structure is largely unknown, without significant additional work. Feb 6 '13 at 23:38

Third answer, which is only a recap of the other answers (this by Aron and this by Christian) present at this time.

For a FEM implementation I would suggest the following steps.

1. Travel element connectivity and compute the sparsity structure of $A$.
2. If a direct solver has to be implemented, renumber d.o.f. to reduce fill-ins, and perform a symbolic factorization of $A$, i.e. compute the sparsity structure of the factors of $A$.
3. Allocate zero sparse matrix in BCRS format with your target sparsity structure (upper triangular part of $A$ for iterative symmetric solver, upper triangular part of Cholesky factor for sparse Cholesky, etc.)
4. Travel elements, compute elemental stiffness matrix, assemble in global sparse matrix $A$.

The important point here is that you have to assemble your matrix only when the sparsity structure is known, so that an efficient storage scheme like BCRS can be used. Of course this is only a hint: you have to choose the sparse matrix storage that better suits your solver.

Wikipedia has an introductory article on sparse matrices, while in Tim David's wonderful site you will fond a treasure of software that can help you in implementing some of these concepts.

Some notes.

Once the sparsity structure is known you can follow Christian's advice for allocating a sparse matrix of zeros in Matlab, and then efficiently assemble your stiffness matrix using standard matlab notation.

The question on how to efficiently compute the sparsity structure of $A$ is still open. I'm not so proficient in Matlab, so my advice is to implement it in python or C, applying some graph theoretical concepts, i.e. consider elements in the matrix as edges in a graph and compute the sparsity structure of the adjacency matrix. List of lists or dictionary of keys are common choices.

• In general you are correct (especially if you need to construct several matrices with the same sparsity structure, e.g., for nonlinear PDEs), but using the sparse constructor as in my answer constructs the sparsity pattern as well as fills the entries at the same time - no need to first build the sparsity pattern and then fill the matrix. (Of course, on the low level, MATLAB is probably doing exactly what you describe.) Feb 7 '13 at 13:27
• @ChristianClason You are right, I've not read your example with the necessary attention. The problem with your approach is memory footprint: you allocate [row, col, val] and then discard it. Typically memory for the sparsity pattern is much less than the memory for the full matrix or for the individual elemental matrices, and can be recycled effectively in the storage of the sparse matrix itself. But, ... inefficient memory allocation is one the most annoying "features" of Matlab itself, so maybe my post is pointless, except for a python or C implementation. Feb 7 '13 at 14:39
• I added my post without reading all the comments to the OP question (they were hidden...) So please forgive if my answer is somehow out of focus. Feb 7 '13 at 15:12
• No worries; your answer adds some important information about what's going on at a lower level. You are right that unnecessary temporary storage is one of the disadvantages of MATLAB; but I just want to point out that row,col,val are arrays and thus only need O(nnz) memory, where nnz is the number of nonzero entries of the sparse matrix. Asymptotically, the memory requirements are therefore not larger than for the sparse matrix itself. The only saving is that the sparsity pattern entries can be stored as bool, which needs less memory than the double for the matrix entries. Feb 7 '13 at 15:24