在 Eigen 中构造稀疏三对角矩阵答案

【问题标题】：Constructing sparse tridiagonal matrix in Eigen在 Eigen 中构造稀疏三对角矩阵
【发布时间】：2018-08-15 02:51:27
【问题描述】：

如何在 Eigen 中构造稀疏的三对角矩阵？我要构造的矩阵在 Python 中如下所示：

alpha = 0.5j/dx**2
off_diag = alpha*np.ones(N-1)
A_fixed = sp.sparse.diags([-off_diag,(1/dt+2*alpha)*np.ones(N),-off_diag],[-1,0,1],format='csc')

如何使用 Eigen 包在 C++ 中完成此操作？看起来我需要使用 here 中记录的“三元组”，但考虑到这应该是一个相当常见的操作，有没有更简单的方法可以做到这一点？

另一个问题是我应该使用行专业还是列专业。我想求解矩阵方程Ax=b，其中A 是一个三对角矩阵。当我们手动进行矩阵向量乘法时，我们通常将矩阵的每一行乘以列向量，因此以行优先存储矩阵似乎更有意义。但是电脑呢？如果要解决Ax=b，首选哪一个？

谢谢

【问题讨论】：

标签： c++ sparse-matrix linear-algebra eigen

【解决方案1】：

三元组是建立稀疏矩阵的指定方法。

您可以采用更直接的方法并使用A.coeffRef(row, col) = val 或A.inser(row,col) = val，即逐个元素填充矩阵。由于您有一个三对角系统，因此您事先知道矩阵的非零数，并且可以使用A.reserve(Nnz) 保留空间。一个愚蠢的方法，但仍然有效，是：

uint N(1000);
CSRMat U(N,N);
U.reserve(N-1);

for(uint j(0); j<N-1; ++j)
    U.insert(j,j+1) = -1;
CSRMat D(N,N);
D.setIdentity();
D *= 2;
CSRMat A = U + CSRMat(U.transpose()) + D;

至于求解器和首选存储顺序，据我回忆，这是次要的。虽然 C(++) 以行优先格式存储连续数据，但是否以最佳方式访问数据取决于算法（行优先存储顺序为逐行）。通常，算法的正确性不取决于数据的存储顺序。它的性能取决于存储顺序和实际数据访问模式的兼容性。如果您打算使用 Eigen 自己的求解器，请坚持使用其默认选择（col-major）。如果您打算与其他库（例如 ARPACK）交互，请选择库喜欢/需要的存储顺序。

【讨论】：