How to assemble large sparse matrices effectively in python/scipy

Question

I am working on an FEM project using Scipy. Now my problem is, that the assembly of the sparse matrices is too slow. I compute the contribution of every element in dense small matrices (one for each element). For the assembly of the global matrices I loop over all small dense matrices and set the matrice entries the following way:

[i,j] = someList[k][l]
Mglobal[i,j] = Mglobal[i,j] + Mlocal[k,l]

Mglobal is a lil_matrice of appropriate size, someList maps the indexing variables.

Of course this is rather slow and consumes most of the matrice assembly time. Is there a better way to assemble a large sparse matrix from many small dense matrices? I tried scipy.weave but it doesn't seem to work with sparse matrices

Answer 1

I posted my response to the scipy mailing list; stack overflow is a bit easier to access so I will post it here as well, albeit a slightly improved version.

The trick is to use the IJV storage format. This is a trio of three arrays where the first one contains row indicies, the second has column indicies, and the third has the values of the matrix at that location. This is the best way to build finite element matricies (or any sparse matrix in my opinion) as access to this format is really fast (just filling an an array).

In scipy this is called coo_matrix; the class takes the three arrays as an argument. It is really only useful for converting to another format (CSR os CSC) for fast linear algebra.

For finite elements, you can estimate the size of the three arrays by something like

size = number_of_elements * number_of_basis_functions**2

so if you have 2D quadratics you would do number_of_elements * 36, for example. This approach is convenient because if you have local matricies you definitely have the global numbers and entry values: exactly what you need for building the three IJV arrays. Scipy is smart enough to throw out zero entries, so overestimating is fine.