在 Eigen 中使用 Cholesky 方法求解大型稀疏线性系统答案

【问题标题】：Solving a large sparse linear system using Cholesky Method in Eigen在 Eigen 中使用 Cholesky 方法求解大型稀疏线性系统
【发布时间】：2019-04-06 09:42:29
【问题描述】：

我实际上是在尝试使用 Eigen 中的 Cholesky 方法解决大型稀疏线性系统。我从this site 下载了一个稀疏矩阵（cfd1）。我试图解决线性系统 Ax = b 定义如下： A 是矩阵 cfd1，b = A * xe 其中 xe 是具有相同大小的矩阵行数的向量。总之，使用 Matlab 符号我想解决：x = A\b。代码如下：

#include <iostream>
#include <Eigen/Dense>
#include <unsupported/Eigen/SparseExtra>
#include<Eigen/SparseCholesky>

using namespace std;
using namespace Eigen;
int main()
{

   SparseMatrix<double> mat;
   VectorXd x;
   loadMarket(mat, "Path of downloaded matrix");

   cout << "Number of Rows:\n" << mat.rows() << endl;

   ArrayXd xe = ArrayXd::Constant(mat.rows(), 1);
   cout << xe << endl;
   SparseVector<double> b = mat*xe;

   SimplicialLLT<SparseMatrix<double> > solver;
   x = solver.compute(mat).solve(b);
   cout << x << endl;



}

问题是当我编译时我得到了当前的错误：

error: invalid operands to binary expression
      ('SparseMatrix<double>' and 'Eigen::ArrayXd' (aka 'Array<double, Dynamic,
      1>'))
   SparseVector<double> b = mat*xe;
                            ~~~^~~
/Users/anto/Desktop/example/eigen-eigen-323c052e1731/Eigen/src/SparseCore/../plugins/CommonCwiseBinaryOps.h:50:29:
note: 
      candidate function template not viable: no known conversion from
      'Eigen::ArrayXd' (aka 'Array<double, Dynamic, 1>') to 'const
      Eigen::SparseMatrixBase<Eigen::SparseMatrix<double, 0, int>
      >::StorageBaseType' (aka 'const
      Eigen::SparseMatrixBase<Eigen::SparseMatrix<double, 0, int> >') for 2nd
      argument

谁能帮我解决它？

【问题讨论】：

标签： c++ sparse-matrix eigen

【解决方案1】：

两件事：

您不能在产品表达式中混合使用Array 和Matrix，即将ArrayXd 替换为VectorXd。
稀疏矩阵与稠密向量的乘积通常是稠密的，因此您需要将乘积分配给稠密的VectorXd 而不是SparseVector

以下编译

SparseMatrix<double> mat;
loadMarket(mat, "Path of downloaded matrix");

cout << "Number of Rows:\n" << mat.rows() << endl;

VectorXd xe = VectorXd::Constant(mat.rows(), 1);
cout << xe << endl;
VectorXd b = mat*xe;

SimplicialLLT<SparseMatrix<double> > solver;
VectorXd x = solver.compute(mat).solve(b);
cout << x << endl;

【讨论】：