Eigen建立稀疏矩阵和迭代求解非对称稀疏矩阵

定义稀疏矩阵：

稀疏矩阵num行num列

SparseMatrix<double> A_sparse(num,num);

稀疏矩阵需要一个储存其元素所在行列位置和元素值的Triplet类型的变量，用std::vector储存起来：

std::vector<Eigen::Triplet<double>> tripletlist;

用push_back把行列信息和对应值填入tripletlist中

for(int i = 0; i != num; ++i) {

        tripletlist.push_back(Eigen::Triplet<double>(i, i, i+1));

        b(i) = i + 1;

    }

    for(int i = 0; i != num; ++i) {

        for(int j = 0; j != num; ++j) {

            if(i + j == num-1)

                tripletlist.push_back(Eigen::Triplet<double>(i, j, i*j+1));

        }

    }

用setFormTriplets将信息填入稀疏矩阵

A_sparse.setFromTriplets(tripletlist.begin(), tripletlist.end());

将稀疏矩阵转化为压缩格式

A_sparse.makeCompressed();

以上就是建立稀疏矩阵A_sparse的流程

下边开始构建求解器，

求解Ax=b，其中x为待求项

求解器类型BiCGSTAB，使用迭代方式求解稀疏矩阵，用于求解方阵。

BiCGSTAB<SparseMatrix<double>> solver;

设置残差

solver.setTolerance(1e-8);

计算矩阵的广义特征值分解

solver.compute(A_sparse);

如果分解失败，这时候

solver.info()!=Success

迭代求解x

x = solver.solve(b);

如果求解失败，这时候

solver.info()!=Success

我们可以这样来计算相对误差

其中的.array()指代对向量的每一个对应元素单独处理

例如x.array() / y.array()等同于x./y，对应元素相除。

VectorXd error = (A_sparse * x - b).array() / b.array();

double maxerror = -10, maxi;

for (int i = 0; i != error.size(); ++i) {

if(fabs(error(i)) > maxerror) {

maxerror = fabs(error(i));

maxi = i;

}

}

完整程序：