C语言实现稀疏矩阵运算，求解线性方程组问题（基于eigen）

近期导师让我将一个matlab代码改写成C，其中涉及到大量的矩阵运算和稀疏矩阵运算，经过一段时间的摸索和实践，记录一下这段时间的学习经验

---

首先由于要基于C语言实现一些矩阵的运算，加减乘除，求逆，以及解线性方程组AX=B之类。

这种矩阵运算自己写的话就有点浪费时间，而且效率肯定不会很高了...故选择eigen库，开源库不用考虑引用问题，且数据结构、函数等相对简单，安装方便，很好上手。

Windows下Eigen + vs配置_wilsonass的博客-CSDN博客

基础的一些矩阵定义及基础矩阵操作可参考以下：

Eigen: Main Page（官方说明文档）

C++矩阵库 Eigen 简介 - rainbow70626 - 博客园 (cnblogs.com)

---

基本求解器（只需要eigen就可以）

下面具体说一下解系数为稀疏矩阵的线性方程组的问题：

在eigen中，有自带的一部分求解器，参考这篇博客的讲解：

Eigen教程5 - 求解稀疏线性方程组_张学志の博客-CSDN博客

如SimplicialLLT、SimplicialLDLT等求解器是可以直接使用的：

SimplicialLDLT<SparseMatrix<float>> solver;
solver.compute(A);
VectorXf x = solver.solve(B);

这是直接法求解，此外还有迭代法求解线性方程组：

87 Eigen应用：解线性方程组的古典迭代法 - 十步一杀2017 - 博客园 (cnblogs.com)

共轭梯度求解

//共轭梯度求解器ConjugateGradient
ConjugateGradient<SparseMatrix<float>, Lower | Upper> cg;
cg.compute(A);
VectorXf x = cg.solve(B);
std::cout << "#iterations:     " << cg.iterations() << std::endl;
std::cout << "estimated error: " << cg.error() << std::endl;
// update b, and solve again
x = cg.solve(B);

最小二乘共轭梯度求解

//最小二乘共轭梯度求解器LeastSquaresConjugateGradient
LeastSquaresConjugateGradient<SparseMatrix<float> > lscg;
lscg.compute(A);
VectorXf x = lscg.solve(B);
std::cout << "#iterations:     " << lscg.iterations() << std::endl;
std::cout << "estimated error: " << lscg.error() << std::endl;
// update b, and solve again
x = lscg.solve(B);

一些求解器的用法参考官方说明文档，里面部分会给出示例代码。

---

速度优化（需要其他库）：

Suitesparse

SuiteSparse是较多用来专门进行稀疏矩阵运算的库，包括矩阵的LU分解，QR分解，Cholesky分解等。其包含了众多的依赖库，例如：blas库、lapack库、cholmod库等，所以安装很复杂。不过值得庆幸的是，国外早有大牛已经实现了在windows,linux或者mac等所平台上的cmake脚本，具体安装过程可参考：

win10+vs2015配置suitesparse1.3.0_玉树银花冬飞雪的博客-CSDN博客

安装好后就可以使用CholmodSimplicialLDLT、CholmodSimplicialLLT等求解器，实测求解速度会快7、8倍。

CholmodSimplicialLDLT<SparseMatrix<double>> solver;
solver.compute(A);
VectorXd x = solver.solve(B);

inter mkl

英特尔的数学核心函数库，是一套经过高度优化和广泛线程化的数学例程，专为需要极致性能的科学、工程及金融等领域的应用而设计。核心数学函数包括 BLAS、LAPACK、ScaLAPACK1、稀疏矩阵解算器、快速傅立叶转换、矢量数学及其它函数。具体安装过程可参考：

vs2015+eigen+intel MKL_pukitoto的博客-CSDN博客_eigen mkl

安装好后就可以使用PardisoLLT，PardisoLDLT，PardisoLU等求解器，在自己测试的求解稀疏矩阵线性方程组的速度上要快于使用SuiteSparse库。

PardisoLDLT<SparseMatrix<double>> solver;
solver.compute(A);
VectorXd x = solver.solve(B);

虽然最后速度以及比一开始缩短很多，但比之matlab还是有1-2倍时间的差距，根据稀疏矩阵的特性应该还有更加快速的求解方法。