问题：求解Ax=b的解，其中，A为大型稀疏矩阵，长和宽分别为1500^2。

一般思路：

（1）A为对称正定矩阵，对A使用cholesky分解。

（2）A为对称不定矩阵，使用LDL‘分解，即：

PAP'=LDL'

其中，L为单位下三角矩阵，D由阶数为1或者2的对角块构成，P是置换矩阵。

（3）不对称矩阵：LU分解。

（4）长方形矩阵（长>宽）：QR分解或者两边同乘以A’，构建对称矩阵，即：

A'Ax=A'b

这里使用LU分解：

流程：

（1）安装MKL。

（2）新建VS工程，设置mkl并行，添加Eigen引用。

（3）头文件定义宏：#define EIGEN_USE_MKL_ALL 表示使用MKL加速。

（4）引入头文件：

#include <Eigen\Sparse> //稀疏矩阵

#include <Eigen\PardisoSupport> //mkl支持

（5）定义并赋值稀疏矩阵A和b。

（6）主函数相关设置：

#include <windows.h>

SYSTEM_INFO sysInfo;

GetSystemInfo(&sysInfo);

std::cerr << "CPU的数目是" << sysInfo.dwNumberOfProcessors << std::endl; //获取cpu的数目

Eigen::initParallel(); //初始化Eigen并行环境

mkl_set_dynamic(sysInfo.dwNumberOfProcessors); //使用mkl把所有的CPU给用上

（6）求解：

Eigen::PardisoLU<Eigen::SparseMatrix<imageDataType>> solver; //定义mkl求解器

solver.analyzePattern(A);

solver.factorize(A);

Eigen::SparseMatrix<imageDataType, Eigen::ColMajor> result = solver.solve(in); //求解结果

本机12核，对于一个(1500^2)*(1500^2)大小的double型矩阵，处理时间9秒。相比mkl加速前，快了3倍，不过领导对于计算时间还是不满意，要上GPU，用流式方法求解。加入后面能实现的话，再贴代码吧……