點點滴滴

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 (減少注入元)，稀疏矩陣的注入元是指執(zhí)行算法后從初始的零值變?yōu)榉橇阒档哪切┰亍?/p>

UMFPACK現(xiàn)在最新的版本是5.6。UMFPACK可以再matlab 中使用，也可以用到C程序中。接下來我們主要介紹如何在C程序中使用。

如何在C中使用UMFPACK

根據(jù)稀疏矩陣的非零元素個數(shù)以及是否實數(shù)，UMFPACK主要提供以下幾個種類的函數(shù)供調用：

1.umfpack_di_*:  非零元素的個數(shù)為 int 型，元素為實數(shù)

2.umfpack_dl_*:  非零元素的個數(shù)為 SuiteSparse_long型，元素為實數(shù)

3.umfpack_zi_*:   非零元素的個數(shù)為 int 型，元素為復數(shù)

4.umfpack_zl_*:  非零元素的個數(shù)為 SuiteSparse_long型，元素為復數(shù)

根據(jù)不同的情況可調用不同的函數(shù)，接下來我們就來看看 上述函數(shù)中的 * 都代表什么。主要有以下5個函數(shù)，也是用來求解Ax=b 的5個步驟：

1.umfpack_*_symbolic:

該函數(shù)會返回一個指向Symbolic類型的void * 指針。主要進行符號分解。

2.umfpack_*_numeric:

該函數(shù)主要對之間所說的矩陣（PAQ，PRAQ，PR-1AQ ）進行數(shù)值分解，需要用到umfpack_*_symbolic返回的結果。

最后返回一個指向Numeric類型的void * 指針。

3.umfpack_*_solve:

該函數(shù)用來求解線性系統(tǒng)（Ax=b，A^Tx=b），會用到umfpack_*_numeric返回的結果。矩陣A必須是方的。

4.umfpack_*_free_symbolic:

釋放umfpack_*_symbolic得到的Symbolic對象。

5.umfpack_*_free_numeric:

釋放umfpack_*_numeric得到的Numeric對象。

在程序中調用上述方法即可，最后我們會給出一個例子。

UMFPACK中稀疏矩陣的表示方法

  現(xiàn)在我們來看一下UMFPACK中稀疏矩陣是如何表示的，UMFPACK中采用compressed column form的格式表示一個稀疏矩陣。

  假設一個M*N的矩陣，有nz個非零元素，umfpack中這樣表示（int 為例）：

     int  Ap[n+1];

     int Ai[nz];

     double Ax[nz];

1.    Ap 是一個列數(shù)+1 的整型數(shù)組，第一個元素為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) ;
}