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倍时间的差距,根据稀疏矩阵的特性应该还有更加快速的求解方法。