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:
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.
Implementing CSC addition and mapping the local stiffness matrix for each element to a empty global matrix then adding successively.
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.