Assemble 来自较小稀疏矩阵的 eigen3 稀疏矩阵

Question

我正在组装耦合多物理系统的雅可比矩阵。雅可比矩阵由每个系统对角线上的块矩阵和耦合的非对角线块组成。我发现最好 assemble 分开分块，然后用投影矩阵对它们求和以获得完整的雅可比矩阵。伪代码（其中 J[i] 是对角线元素，C[ij] 耦合，P 是对完整矩阵的投影）。

// diagonal blocks
J.setZero();
for(int i=0;i<N;++i){
    J+=P[i]J[i]P[i].transpose()
}
// off diagonal elements
for(int i=0;i<N;++i){
    for(int j=i+1;j<N;++j){
        J+=P[i]C[ij]P[j].transpose()
        J+=P[j]C[ji]P[i].transpose()
    }
}

这需要很多性能，大约占整个程序的 20%，这对于某些汇编来说太多了。由于系统是高度非线性的，因此我必须在每个时间步重新计算雅可比矩阵。 Valgrind 表明资源消耗方法是 Eigen::internal::assign_sparse_to_sparse 并且在这个方法中调用 Eigen::SparseMatrix<>::InsertBackByOuterInner.

有没有更有效的方法assemble这样的矩阵？

（我还必须使用 P*(JP.transpose()) 而不是 PJ*J.transpose() 来编译程序, 可能已经有问题了)

P.S: NDEBUG 和优化开启

编辑：通过将 P.transpose 存储在一个额外的矩阵中，我得到了更好的性能，但总和仍然占程序的 15%

Answer 1

通过就地工作，您的代码会更快。首先，估计最终矩阵中每列的非零数并保留 space（如果尚未完成）：

int nnz_per_col = ...;
J.reserve(VectorXi::Constant(n_cols, nnz_per_col));

如果每列的 nnz 数量非常不均匀，那么您也可以按列计算它：

VectorXi nnz_per_col(n_cols);
for each j
  nnz_per_col(j) = ...;
J.reserve(nnz_per_col);

然后手动插入元素：

for each block B[k]
  for each elements i,j
    J.coeffRef(foo(i),foo(j)) += B[k](i,j)

其中 foo 实施适当的索引映射。

并且下一次迭代不需要保留，但是需要在保留结构的同时将系数值设置为零：

J.coeffs().setZero();

Assemble 来自较小稀疏矩阵的 eigen3 稀疏矩阵

Assemble eigen3 sparsematrix from smaller sparsematrices

sparse-matrix

eigen