How to efficiently assemble global stiffness matrix in sparse storage format (c++)

Question

I am writing a finite element solver in C++. The main bottle neck is assembling the global stiffness matrix in sparse compressed row storage (so far I am only solving steady problems). Because I don't know how many nonzero entries exist in each row, I am currently assuming a constant upper bound on the number of non-zeros per row. My codes currently looks like:

int nrow[fem.Np+1];          //row vector in CRS format initialized to 0
int ncol[fem.Np*MAX_NNZ];    //column vector in CRS format initialized to -1
double A[fem.Np*MAX_NNZ];    //values vector in CRS format initialized to 0.0

//update nrow, ncol, and A 
//MAX_NNZ is the constant upper bound on the number of nonzeros per row
//Np=='number of points (nodes)
//Npe=='number of points per element e.g. 3 for linear triangles
//kmat[Npe][Npe] is the previously computed local element stiffness matrix
for(int i=0;i<e.Npe;i++)
{
  for(int j=0;j<e.Npe;j++)
  {
    r = e.nodes[i];
    c = e.nodes[j];
    for(int p=MAX_NNZ*r-MAX_NNZ;p<MAX_NNZ*r;p++)
    {
      if(ncol[p]==-1)
      {
        for(int l=r;l<fem.Np+1;l++)
          nrow[l]++;
        ncol[p] = c-1;
        A[p] = e.kmat[i][j];
        break;
      }
      else if(ncol[p]==c-1)
      {
        A[p] = A[p] + e.kmat[i][j];
        break;
      }
    }
  }

Is this the most efficient way to assemble the global matrix in CRS format? Are there improvements that I can make? Thanks.

Answer 1

When I do assembly, I use a Coordinate sparse matrix format. This is basically just a list of (row, col, value) tuples. This is especially useful for sparse matrix construction when the exact pattern of non-zeros is unknown (or simply not yet computed).

In my implementation, I prohibit myself from actually querying the value of A(i,j), but this allows me to have a constant-time "add" function since the new value is simply appended to the end of the list of tuples. After construction is complete, I sort the list of tuples (summing any duplicate entries). When I want to use matrix-vector operations, I convert from this sorted coordinate-format sparse matrix to CRS. Performing the sort and conversion takes O(N logN) time, bottlenecked by the sorting step.

A different approach to this problem is to use your mesh to pre-compute where the non-zero entries are going to be prior to starting assembly. This allows you to efficiently update the values of the nonzeros. Trying to assemble directly in CRS without doing this step first is generally considered a bad idea, since you will have to move a lot of data around if you introduce new nonzeros.

Answer 2

Just as a reference point, in the deal.II finite element library (see http://www.dealii.org), we typically pre-compute the sparsity pattern once at the beginning by "simulating" what would happen if we did the assembly, and then create the sparse matrix from this sparsity pattern.

The first tutorial in which we build a sparsity pattern is step-3: http://www.dealii.org/developer/doxygen/deal.II/step_3.html We do this in the Step3::setup_system() function. Sparsity patterns in general are discussed here: http://www.dealii.org/developer/doxygen/deal.II/group__Sparsity.html

Answer 3

Have a look at the thread below. Accepted answer gives a good explanation of how to create a CRS structure on the fly.

https://stackoverflow.com/questions/13026111/fast-accessing-elements-of-compressed-sparse-row-csr-sparse-matrix

Answer 4

To do that, I am using two different sparse matrix data structures. One is dynamic and new coefficients can be inserted into it during assembly. Once assembly is finished, I convert it into a CRS matrix. It is not as optimal as precomputing the sparsity pattern with two passes, but easier to use, and works reasonably well (I used it for problems of dimension 20 million with hundreds of millions non-zero entries).

Implementation is available in my OpenNL and Geogram libraries:

http://alice.loria.fr/software/geogram/doc/html/index.html

http://alice.loria.fr/index.php/software/4-library/23-opennl.html

My dynamic matrix structure is an array of SparseRow objects. The idea is to have objects that have both a size (number of used coefficients) and a capacity (number of allocated coefficients), and to double capacity each time more space is needed (to avoid reallocating too often). Same strategy is used in the std::vector C++ class:

/**
 * \brief Represents a coefficient in a sparse matrix
 * \relates NLSparseMatrix
 */
typedef struct  {
    /**
     * \brief index of the coefficient.
     */    
    NLuint index ;

    /**
     * \brief value of the coefficient. 
     */    
    NLdouble value ; 
} NLCoeff ;

/**
 * \brief Represents a row or a column of a sparse matrix
 * \relates NLSparseMatrix
 */
typedef struct {
    /**
     * \brief number of coefficients. 
     */    
    NLuint size ;

    /** 
     * \brief number of coefficients that can be 
     * stored without reallocating memory.
     */    
    NLuint capacity ;

    /**
     * \brief the array of coefficients, with enough
     * space to store capacity coefficients.
     */
    NLCoeff* coeff ;  
} NLRow ;

There are functions to construct a row, to add a new coefficient to it and to delete it:

void nlRowConstruct(NLRow* c) {
    c->size     = 0 ;
    c->capacity = 0 ;
    c->coeff    = NULL ;
}

void nlRowDestroy(NLRow* c) {
    free(c->coeff) ;
}

void nlRowGrow(NLRow* c) {
    if(c->capacity != 0) {
        c->capacity = 2 * c->capacity ;
        c->coeff = (NLCoeff*)realloc(c->coeff, c->capacity*sizeof(NLCoeff)) ;
    } else {
        c->capacity = 4 ;
        c->coeff = (NLCoeff*)malloc(c->capacity*sizeof(NLCoeff)) ;
    }
}

void nlRowAdd(NLRow* c, NLuint index, NLdouble value) {
    NLuint i ;
    for(i=0; i<c->size; i++) {
        if(c->coeff[i].index == index) {
            c->coeff[i].value += value ;
            return ;
        }
    }
    if(c->size == c->capacity) {
        nlRowColumnGrow(c) ;
    }
    c->coeff[c->size].index = index ;
    c->coeff[c->size].value = value ;
    c->size++ ;
}

A SparseMatrix is like:

typedef struct {
   /**
     * \brief number of rows 
     */    
    NLuint m ;

    /**
     * \brief number of columns 
     */    
    NLuint n ;

    /**
      * \brief the rows, array of size m
      */
        NLRow* rows;
    } SparseMatrix;

Inserting a coefficient in the matrix is done as follows:

void nlSparseMatrixAdd(NLSparseMatrix* M, NLuint i, NLuint j, NLdouble value) {
    nl_parano_range_assert(i, 0, M->m - 1) ;
    nl_parano_range_assert(j, 0, M->n - 1) ; 
    nlRowColumnAdd(&(M->row[i]), j, value) ;  
}