# Fast way to build stiffness directly as CSC matrix

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.

• I am pretty sure that you need to use HB or RB matrix formats (which are CSC type formats) in elemental mode but I cannot find the reference I am looking for. So I am going to give you closest I can. people.math.sc.edu/Burkardt/data/rb/rb.html . This is the format I used to move files around when I was doing my master's research – Abdullah Ali Sivas Jun 3 '20 at 14:28
• Also you might want to check the documentation for PETSc AIJ matrices and similarly hypre matrices. I think they use a COO-like format but I honestly don't remember. It has been at least four years since I last wrote a FEM assembly code in PETSc. – Abdullah Ali Sivas Jun 3 '20 at 14:34
• What does COO stand for? – Wolfgang Bangerth Jun 3 '20 at 16:36
• According to SciPy sparse library, COO, DOK and LIL are recommended for this purpose. I have used COO for FEM Py thon and it seems to work, but maybe I have not used large enough problems. – nicoguaro Jun 3 '20 at 19:43

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})$$.

• oh wow I should have looked into your source code. I have step-35 running as we speak! This is also a great answer. Will try this now. I understand now that the sparsity pattern for DOFhandler is actually a by-product of your stiffness matrix build process. Vielen Dank! – Not a chance Jun 3 '20 at 19:49
• @Notachance The (calling side of this) code is here: dealii.org/developer/doxygen/deal.II/… Everything in blue is clickable. You should find all three steps from above represented. – Wolfgang Bangerth Jun 3 '20 at 22:10
• This was helpful. My code is much faster. The limiting factor on assembly time is now integration rather than juggling sparse matrix entries. Thanks again. – Not a chance Jun 5 '20 at 14:15
• @Notachance -- fantastic, that's what should happen! – Wolfgang Bangerth Jun 5 '20 at 20:22