稀疏矩阵解题数学库 -- UMFPACK

UMFPACK就是求解类似于Ax=b这样问题的一个库，来自佛罗里达州立大学。可以直接到http://www.cise.ufl.edu/research/sparse/umfpack/去下载对应的包然后编译得到Windows下的lib，包含到自己的工程就可以了。 

什么是UMFPACK

UMFPACK是专门求解类似于Ax=b这样的稀疏矩阵方程的一个库，一般情况下A 是稀疏非对称的矩阵。它是基于非对称的多波前算法（ Unsymmetric-pattern MultiFrontal method ）对稀疏矩阵方程进行求解的。UMFPACk可以对PAQ，PRAQ，PR-1AQ 这样的矩阵进行 LU分解（L 和 U 分别是下三角矩阵和上三角矩阵，P和Q是可置换矩阵，R是对角阵） 。这里有一个概念：reduce fill-in (减少注入元)，稀疏矩阵的注入元是指执行算法后从初始的零值变为非零值的那些元素。

UMFPACK现在最新的版本是5.6。UMFPACK可以再matlab 中使用，也可以用到C程序中。接下来我们主要介绍如何在C程序中使用。

如何在C中使用UMFPACK

根据稀疏矩阵的非零元素个数以及是否实数，UMFPACK主要提供以下几个种类的函数供调用：

1.umfpack_di_*:  非零元素的个数为 int 型，元素为实数

2.umfpack_dl_*:  非零元素的个数为 SuiteSparse_long型，元素为实数

3.umfpack_zi_*:   非零元素的个数为 int 型，元素为复数

4.umfpack_zl_*:  非零元素的个数为 SuiteSparse_long型，元素为复数

根据不同的情况可调用不同的函数，接下来我们就来看看 上述函数中的 * 都代表什么。主要有以下5个函数，也是用来求解Ax=b 的5个步骤：

1.umfpack_*_symbolic:

该函数会返回一个指向Symbolic类型的void * 指针。主要进行符号分解。

2.umfpack_*_numeric:

该函数主要对之间所说的矩阵（PAQ，PRAQ，PR-1AQ ）进行数值分解，需要用到umfpack_*_symbolic返回的结果。

最后返回一个指向Numeric类型的void * 指针。

3.umfpack_*_solve:

该函数用来求解线性系统（Ax=b，A^Tx=b），会用到umfpack_*_numeric返回的结果。矩阵A必须是方的。

4.umfpack_*_free_symbolic:

释放umfpack_*_symbolic得到的Symbolic对象。

5.umfpack_*_free_numeric:

释放umfpack_*_numeric得到的Numeric对象。

在程序中调用上述方法即可，最后我们会给出一个例子。

UMFPACK中稀疏矩阵的表示方法

  现在我们来看一下UMFPACK中稀疏矩阵是如何表示的，UMFPACK中采用compressed column form的格式表示一个稀疏矩阵。

  假设一个M*N的矩阵，有nz个非零元素，umfpack中这样表示（int 为例）：

     int  Ap[n+1];

     int Ai[nz];

     double Ax[nz];

1.    Ap 是一个列数+1 的整型数组，第一个元素为0 ，即Ap[0]=0; 之后的每个元素为当前列与之前所有列的非零元素之和。
2.    Ai 对应每列中非零元素所在的位置，Ai 依赖于Ap。
3.    Ax 与Ai 对应，表示每个非零元素的值。 

如上所示一个稀疏矩阵，对应的Ap, Ai, Ax如下所示:

 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};

小例子demo

#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) ; }