定義稀疏矩陣：

?

稀疏矩陣num行num列

SparseMatrix<double> A_sparse(num,num);

?

稀疏矩陣需要一個儲存其元素所在行列位置和元素值的Triplet類型的變量，用std::vector儲存起來：

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

?

用push_back把行列信息和對應(yīng)值填入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());

?

將稀疏矩陣轉(zhuǎn)化為壓縮格式

A_sparse.makeCompressed();

?

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

?

下邊開始構(gòu)建求解器，

求解Ax=b，其中x為待求項

?

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

BiCGSTAB<SparseMatrix<double>> solver;

?

設(shè)置殘差

solver.setTolerance(1e-8);

?

計算矩陣的廣義特征值分解

solver.compute(A_sparse);

?

如果分解失敗，這時候

solver.info()!=Success

?

迭代求解x

x?=?solver.solve(b);

?

如果求解失敗，這時候

solver.info()!=Success

?

我們可以這樣來計算相對誤差

其中的.array()指代對向量的每一個對應(yīng)元素單獨處理

例如x.array() / y.array()等同于x./y，對應(yīng)元素相除。

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;

}

}

?

完整程序：