點點滴滴

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