Efficient assembly of finite element matrix in MATLAB

Question

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.

Answer 1

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.

Answer 2

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

Answer 3

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.

Answer 4

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.