umfpack求解求解器求解大型稀疏矩阵方程详解

计算通常会引入矩阵并得到一个如下的矩阵方程：
Ax = b
其中的x就是一个关于目标点的参数列向量，b是对应的观测值向量，而A就是根据应该背景得到的数据矩阵。
可能在大多数情况下得到的A会是一个规模较大且较为稀疏的矩阵。UMFPACK就是求解类似于Ax=b这样问题的一个库

一个最简单的使用程序如下：

#include "umfpack.h"
int n = 5 ;
int Ap [ ] = {0, 2, 5, 9, 10, 12} ;
int Ai [ ] = { 0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4} ;
double Ax [ ] = {2.0, 3.0, 3.0, -1.0, 4.0, 4.0, -3.0, 1.0, 2.0, 2.0, 6.0, 1.0} ;
double b [ ] = {8.0, 45.0, -3.0, 3.0, 19.0} ;
double x [5] ;
int main()
{
	double *null = (double *) NULL ;
	int i ;
	void *Symbolic, *Numeric ;
	(void) umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null) ;
	(void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
	umfpack_di_free_symbolic (&Symbolic) ;
	(void) umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null) ;
	umfpack_di_free_numeric (&Numeric) ;
	// ... 从x中得到最终的解并使用即可
	return (0) ;
}

其中对应的A原型为：

UMFPACK中求解各种Symboli以及Numeric就不过多关注了，主要来看一下UMFPACK中矩阵的表示方法。一般情况下，对于大型的稀疏矩阵，采用压缩的方法来表示更能节省空间。

什么意思呢？
这里就涉及到一个 稀疏矩阵和矩阵压缩的概念及实现*

就是说把一个矩阵的非零元素取出来形成一个矩阵再进行计算。

详解可以看
https://blog.csdn.net/weixin_42000717/article/details/124704649?spm=1001.2014.3001.5501

UMFPACK里用的压缩方法也比较简单直白，UMFPACK采用列主的矩阵索引：
Ap为一个整形数组，其大小为：A的列数+1，第一个元素为0，之后的每个元素取值为当前列与之前所有列中的非0原素个数之和。如Ap[1]=2表示第一列中非0元素数为2个；Ap[2]=5表示第二列中非0元素个数为5-2=3个；依次类推。
Ai在Ap的背景下使用，对应每列中非0元素所在的位置，比如Ai[0],Ai[1]就表示第一列中的两个非0元素的位置为0,1。
Ax就与Ai一一对应了，表示每个非0元素的具体值。
b与x就是对应的列向量了，常规存储方式。

另外，注意函数umfpack_di_solve中使用的宏UMFPACK_A，它用来说明要求解的矩阵方程的类型，对应关系为：
UMFPACK_A: Ax = b
UMFPACK_At: A’x=b
UMFPAC_Aat: A.'x=b

实际应用中解一个矩阵A的规模为4611x4611的方程，平均每行有七个左右的非0元素（最多9个；最少4个），UMFPACK对该方程的计算时间为0.025ms。