从几个较小的稀疏矩阵构建大型稀疏矩阵答案

【问题标题】：Build large sparse matrix from several smaller sparse matrices从几个较小的稀疏矩阵构建大型稀疏矩阵
【发布时间】：2021-09-22 00:53:46
【问题描述】：

我正在使用 Eigen 库来解决一个 FEM 问题，在该问题中，我使用几个类似的稀疏矩阵来计算我计算的不同类型的导数。为了构建稀疏矩阵来解决系统问题，我想使用逗号初始化器，但稀疏矩阵（https://eigen.tuxfamily.org/dox/group__TutorialSparse.html 和 https://stackoverflow.com/a/43532618/15811117）不支持。在那个答案中，Henri Menke 确实建议使用 Eigen 不受支持的 Kronecker 产品，尽管我认为这不会在这里起作用。

使用 NxM 网格，我计算了维度为 N*M x N*M 的几个导数矩阵（稀疏的）Dx2, Dxz, 和 Dz2。我的求解器矩阵如下，（K 只是常数，0 只是填充物——一个零矩阵，N*M x N*M）

K3*Dx2 + K1*Dz2   K2*Dxz            0              0                 0
K2*Dxz            K3*Dz2 + K1*Dx2   0              0                 0
0                 0                 K1*(Dx2+Dz2)   0                 0
0                 0                 0              K3*Dx2 + K1*Dz2   K2*Dxz
0                 0                 0              K2*Dxz            K3*Dz2 + K1*Dx2

到目前为止，我通过使用逗号初始化器将稀疏转换为密集，然后通过 .sparseView() 转换回稀疏来构建它。唯一的问题是，如果网格远大于 16x16，我会收到 std::bad_alloc 错误（这是合理预期的，因为它创建了一个巨大的矩阵：5*N*M x 5*N*M），我需要一个至少 100x100 的网格（并且不会总是像当前显示的那样多 0 矩阵，尽管它仍然是一个稀疏矩阵）。

构建此矩阵的最佳方法是什么？我曾想过尝试使用 Triplets 来做到这一点，类似于此 (Convert an Eigen matrix to Triplet form C++)。

编辑：我不认为使用 Triplets 会那么容易（尽管我目前正在尝试），因为我在每个 D_ _ 矩阵中为边界条件编辑了几个值......除非有办法在给定行号和列号的 std::vector 中找到 Triplet？

更新： Triplet 和 Kronecker product 方法都是可行的，但我遇到了几个问题：

使用 Kronecker 产品，我在编译时收到以下错误（大约一百次左右）：

.../eigen/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h:225:97: error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<int, std::complex<double>, Eigen::internal::scalar_product_op<int, std::complex<double> > >’

使用 Triplet 方法，我定义了一些运算符和一个函数来将矩阵“移动”到正确的位置（基本上是 Kronecker 产品正在做的事情）。我只需要找出索引越界问题。

【问题讨论】：

导数矩阵是否稀疏？
@CJR 是的，它们很稀疏；它们最初是使用 setFromTriplets 函数制作的。我现在正在尝试看看是否可以将它们保留为三元组的向量，直到我建立这个矩阵。
您为什么认为 Kronecker 产品不是一个可行的解决方案？使用克罗内克积，您可以创建重复 Dx2、Dxz 和 Dz2 获得的矩阵，然后将它们求和以获得最终矩阵。
@LGrementieri 我想我误解了如何使用 Kronecker 产品——我同意，它应该可以工作。可悲的是，我在编译时遇到了数百个错误:( 请参阅更新。

标签： c++ sparse-matrix eigen bad-alloc

【解决方案1】：

问题在于我在 Kronecker 产品中使用的矩阵类型为 SparseMatrix<int> 和 SparseMatrix<complex<double>>。 int 和 complex<double> 没有运算符重载，因此通过使两个矩阵都兼容类型，代码可以按预期编译和运行。

【讨论】：