Fast way to build stiffness directly as CSC matrix

Question

I have been working on finite element code in Fortran 2008, and have implemented my own sparse matrix types. I have found that mapping local stiffness matrices (real type) to a global COO sparse type and then converting to CSC works well but above a certain size, sorting the COO becomes prohibitive.
Instead I would like to build the CSC matrix directly, but I have not been able to make it go fast enough. I don't want to bore anyone with specifics of my code, but I have tried the following:

1. Inserting values from each local stiffness matrix into a predeclared global CSC matrix. This was slow because allocating new space and shifting the column pointers to the right is slow.

2. Implementing CSC addition and mapping the local stiffness matrix for each element to a empty global matrix then adding successively.

3. Tried building an adjacency matrix so that all the entries would be preallocated. This proved to be just as slow.

I'm interested in your structural ideas here. How do people normally do this? My code is in 3D using 2nd Order elements, on an unstructured grid. A smaller problem has ~4 million nodes.

FYI the problem with COO-> CSC transfer on large matrices is that I implemented mergesort, and my implementation crawls once the total data is greater than around 32MB.

Answer 1

COO is an unsuitable matrix format except for particular purposes (e.g., if there is a substantial number of rows that have no entries at all, possibly with the exception of the diagonal).

The way all large-scale codes I know of build the system matrix is through a three-step process:

• In the first step, you loop over all cells and figure out which entries you would write into, and from that build the sparsity pattern of the CSC/CSR matrix. To do this, one generally keeps a sorted set of indices for each row of the sparsity pattern (i.e., not a set of index pairs for the entire matrix).

• You then count for each row how many unique column indices are to be written to, and can build a static (unchangeable) CSC/CSR sparsity pattern in which you keep two arrays: One very long array of column indices, and one array of size equal to the number of rows that tells you where in the long array of column indices, the $i$th row starts.

• With the sparsity pattern built, you then allocate the memory for the matrix. This is an array of length equal to the number of column indices.

In the first step, the sort operations only ever sort column indices in an individual row. Furthermore, for each degree on one cell, you know which columns you want to write into; so sort these column indices and then just merge this sorted array with the sorted array of previously added indices for this row. The merge can be done in $O(n_\text{nonzeros per row})$, whereas the initial sort takes $O(n_\text{dofs per cell} \log n_\text{dofs per cell})$.