如何高效组装 FEM 稀疏矩阵

Question

全部成交， 感谢您花时间阅读我的问题。 我正在使用 Eigen3.3.4(http://eigen.tuxfamily.org/index.php?title=Main_Page) 编写一些 FEM 代码。

我阅读了Eigen3.3.4文档，在这个网站（http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html）中，它说 我们应该使用

Ref<MatrixBase>

来避免额外的副本并获得高性能。

所以在我的 FEM 代码中，对于稀疏矩阵组装部分，假设函数是：

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V,
Ref<SparseMatrix<double> > AMATRIX,Ref<VectorXd> RHS)

其中U代表位移，V代表速度项。 AMATRIX 是我的稀疏矩阵，RHS 是残差项。

然后我尝试在组装之前首先初始化我的 AMATRIX（我有一个三元组列表，其中包含所有非零元素及其值（我将值设置为零以进行初始化）） 所以我尝试了：

AMATRIX.setFromTriplets(ZeroTripList.begin(),ZeroTripList.end());

但是我有一个错误：

class Eigen::Ref<Eigen::SparseMatrix<double, 0, int> >’ has no member named ‘setFromTriplets

那么我该如何解决这个问题呢？

我的解决方案之一是使用：

FormFE(const Ref<VectorXd> &U,const Ref<VectorXd> &V,
SparseMatrix<double> &AMATRIX,Ref<VectorXd> RHS)

这工作得很好，但我不确定它是否有效。 我不太擅长 cpp :P

其实我的问题是：

1. 如何高效使用Eigen（尤其是FEM计算），我在每个FEM相关函数中几乎到处都使用Eigen的VectorXd和MatrixXd。
2. 如何高效组装SparseMatrix？
3. 是否可以对 FEM 组装进行一些 OpenMP 并行化？
4. 欢迎对基于 C++ 的 FEM 编码提出任何有用的建议（库推荐或任何有用的想法）！

谢谢你。 致以诚挚的问候。

Answer 1

是的，通过

SparseMatrix<double> &

是正确的做法。

Ref<SparseMatrix>

的目的是传递类似于

SparseMatrix

的组装对象，例如子稀疏矩阵、

Map<SparseMatrix>

...

使用 setFromTriplets 也是确保获得良好性能的正确做法。如果操作正确（即正确调用保留和正确的插入顺序），使用

mat.insert(i,j) = val;

直接插入元素可能会快 2 倍。但如果你弄错了，它也可能慢 100 倍......请参阅文档。使用

SparseMatrix::insert

也可以使用 OpenMP 填充矩阵，但这需要更加小心和严格，这是一个典型的模式：

int n_cols = ??, n_rows = ??;
std::vector<int> nnz_per_col(n_cols);
// set each nnz_per_col[j] to the exact number
// of non-zero entries in the j-th column (or more, but NOT less)
SparseMatrix<double> mat(n_rows, n_cols);
#pragma omp parallel for
for(int j=0; j<cols; ++j) {
  for each non zero entry i in the j-th column {
    // preferably with increasing i
    double val_i_j = ...;
    mat.insert(i,j) = val_i_j;
  }
}

当然，如果您更容易的话，您也可以按行工作。在这种情况下，请使用

SparseMatrix<double,RowMajor>

。当然，您可以调整此模式以处理列/行块等。

如果对于装配，您需要使用一些密集的矩阵/向量，那么我相信它们非常小且大小固定。那么最好使用静态分配的

Matrix<double,N,M>

和

Matrix<double,N,1>

类型，而不是使用 MatrixXd/VectorXd。这将防止大量内存分配/释放。

最后，最重要的建议：如果您关心性能，请不要忘记在研究优化代码所需的时间和精力之前分析您的代码。另外，始终在编译器优化打开的情况下进行基准测试/分析。

Answer 2

对于同样的问题，我并行创建了 std::vectorEigen::Triplet （每个线程并行收集到线程本地向量，最后我合并这些向量）。并行组装三元组的最有效方法是什么？我是否必须对每列行进行排序，合并相同的出现，然后将每列分配给一个线程？有没有办法避免排序？