当我进行组装时，我使用坐标稀疏矩阵格式。这基本上只是一个 (row, col, value) 元组的列表。当非零的确切模式未知（或根本尚未计算）时，这对于稀疏矩阵构造特别有用。

在我的实现中，我禁止自己实际查询 A(i,j) 的值，但这允许我拥有一个恒定时间的“add”函数，因为新值只是简单地附加到元组列表的末尾。构建完成后，我对元组列表进行排序（对所有重复条目求和）。当我想使用矩阵向量运算时，我从这个排序的坐标格式稀疏矩阵转换为 CRS。执行排序和转换需要 O(N logN) 时间，这是排序步骤的瓶颈。

解决此问题的另一种方法是使用您的网格在开始组装之前预先计算非零条目的位置。这使您可以有效地更新非零值。尝试直接在 CRS 中组装而不先执行此步骤通常被认为是一个坏主意，因为如果您引入新的非零值，您将不得不移动大量数据。

作为参考，在 deal.II 有限元库（参见http://www.dealii.org）中，我们通常在开始时通过“模拟”如果我们执行组装，然后从这个稀疏模式创建稀疏矩阵。

我们构建稀疏模式的第一个教程是第 3 步：http: //www.dealii.org/developer/doxygen/deal.II/step_3.html 我们在 Step3::setup_system() 函数中执行此操作。一般的稀疏模式在这里讨论：http: //www.dealii.org/developer/doxygen/deal.II/group__Sparsity.html

看看下面的线程。接受的答案很好地解释了如何动态创建 CRS 结构。

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

为此，我使用了两种不同的稀疏矩阵数据结构。一种是动态的，可以在组装过程中插入新的系数。组装完成后，我将其转换为 CRS 矩阵。它不如预先计算两次通过的稀疏模式最佳，但更易于使用，并且工作得相当好（我将它用于具有数亿非零条目的 2000 万维问题）。

我的 OpenNL 和 Geogram 库中提供了实现：

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

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

我的动态矩阵结构是一个 SparseRow 对象数组。这个想法是让对象同时具有大小（使用的系数的数量）和容量（分配的系数的数量），并且每次需要更多空间时将容量加倍（以避免过于频繁地重新分配）。在 std::vector C++ 类中使用了相同的策略：

/**
 * \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 ;

有一些函数可以构造一行，向它添加一个新系数并删除它：

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++ ;
}

一个稀疏矩阵就像：

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

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

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

在矩阵中插入一个系数如下：

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) ;  
}