UMFPACK调用的接口

在调用UMFPACK的过程中，只需要关心Ap Ai Ax的产生，实现其过程分为下面两种方法：

（1）通过Eigen库，先让矩阵A以稀疏矩阵格式存储（知道矩阵A的非零元素的分布）

//#include <Eigen/Eigen>
#include <Eigen/Sparse>
//#include <umfpack.h>
//#include <Eigen/src/UmfPackSupport/UmfPackSupport.h>
//注意：只有debug版本调试
//参考资料：科学计算中的偏微分方程有限差分法 张文生 高等教育出版社
//                4.7节 边界条件的处理 4.7.2 Neumann边界 P155
// 主要是形成Eigen中要求的稀疏矩阵 而非一般的二维数组
#include <iostream>
#include <vector>

using namespace Eigen;

using namespace std;

typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double

int main()
{
        
    /*---以下几行作为测试用---*/
    Eigen::Vector2d v1, v2;     //Eigen中的变量
    v1 << 5, 6;   //默认的向量为列向量
    cout  << "v1 = " << endl << v1 << endl;
    v2 << 4, 5 ;
    Matrix2d result = v1*v2.transpose();
    cout << "result: " << endl << result << endl;
    cout<<"test----"<<7%int(3.0)<<endl;
    cout<<"Please input the dimension of the Matrix----( N)!"<<endl;
    int N=9;

    SpMat A(N,N); 
    typedef Eigen::Triplet<double> Tri;
    vector<Tri> coefficients;
    //coefficients.push_back(Tri(0,0,2.0));
    for(int i=0;i<sqrt( double(N) );i++)
        coefficients.push_back(Tri(i,i,1.0));
    for(int i=int(sqrt( double(N) ));i<N-sqrt( double(N) );i++)
    {
        if(i%int(sqrt( double(N) ))==0) //B矩阵中的首行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i+1,-2.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-1.0));
            coefficients.push_back(Tri(i,i+int(sqrt(double(N))),-1.0));
        }
        else if((i-(int(sqrt(double(N)))-1))%int(sqrt( double(N) ))==0) //B矩阵中的末行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i-1,-2.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-1.0));
            coefficients.push_back(Tri(i,i+int(sqrt(double(N))),-1.0));
        }
        else //B矩阵的中间行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i-1,-1.0));
            coefficients.push_back(Tri(i,i+1,-1.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-1.0));
            coefficients.push_back(Tri(i,i+int(sqrt(double(N))),-1.0));
        }
    }
    for(int i=N-int(sqrt( double(N) ));i<N;i++)
    {
        if(i==N-int(sqrt( double(N) ))) //B矩阵中的首行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i+1,-2.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-2.0));
        }
        else if(i==N-1) //B矩阵中的末行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i-1,-2.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-2.0));
        }
        else //B矩阵的中间行
        {
            coefficients.push_back(Tri(i,i,4.0));
            coefficients.push_back(Tri(i,i-1,-1.0));
            coefficients.push_back(Tri(i,i+1,-1.0));
            coefficients.push_back(Tri(i,i-int(sqrt(double(N))),-2.0));
        }
    }

    A.setFromTriplets(coefficients.begin(),coefficients.end());
    cout<<endl;
    int _index=A.nonZeros();
    cout<<_index<<endl;
    
    int n=A.cols();
    int *Ap=new int[n+1];
    Ap[0]=0;
    int num=A.nonZeros();
    int *Ai=new int[num];
    double *Ax=new double[num];

    int k=0;
    for(int i=0;i<A.outerSize();i++)
    {
        Ap[i+1]=Ap[i];
        for (Eigen::SparseMatrix<double>::InnerIterator it(A,i); it; ++it)
        {
            Ax[k]=it.value();
            
            Ai[k]=it.row();
            //cout<<Ai[k]<<" ";
            k++;
            Ap[i+1]++;
        }
        cout<<Ap[i]<<" ";
    }
cout<<Ap[A.outerSize()]<<" ";

    return 0;
}

 

（2）矩阵A为二维数组

//Data:2013-2-24
//修改了Ai Ax的类型 利用最大维数n*n来保存，可以调用正确结果 不过不知系统随机分配的值 函数没有用
//Data:2013-2-26
//调用UMFPACK包来实现求解方程组
//UMFPACK采用CSC（列压缩存储） matlab中的接口为A/b
#include <stdio.h>
#include <math.h>
#include "umfpack.h"
//传递的四个参数A b x n -----Data:2013-02-27
//意思为Ax=b n为矩阵维数
int umf(double **A,double *b,double *x,int n)
{
    //printf("--");
    //printf("%d \n",A[1][1]);
    // int n=5;
    double *null =(double *)NULL ;
    void *Symbolic, *Numeric ;
    int i,j;

/*    
    //定义矩阵A
    int **A=new int *[n];
    for(j=0;j<n;j++)
        A[j]= new int[n];

    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            scanf("%d",&A[i][j]);
*/

    int *Ap=new int [n+1];
    Ap[0]=0;

   /* int *Ai=new int[n*n];
    for(i=0;i<n*n;i++)
        Ai[i]=1;
    double *Ax=new double[n*n];
    for(i=0;i<n*n;i++)
        Ax[i]=1;*/
    double epsilon=0.00001;
    int NZnum=0;//矩阵非零元的个数

    int k=0,l=0,m=0;
    for(j=0;j<n;j++)
    {
        for(i=0;i<n;i++)
        {
            if(abs(A[i][j])>epsilon)
            {
             /*   Ai[l++]=i;
                Ax[m++]=A[i][j];*/


                NZnum++;
            }
        }
        Ap[k+1]=NZnum;
        k++;
    }
    Ap[n]=NZnum;

     int *Ai=new int[NZnum];
    double *Ax=new double[NZnum];

    // int k=0,l=0,m=0;
    for(j=0;j<n;j++)
    {
        for(i=0;i<n;i++)
        {
            if(abs(A[i][j])>epsilon)
            {
                Ai[l++]=i;
                Ax[m++]=A[i][j];
            }
        }
    }

    //printf("--Ai---\n");
    //for (int i=0; i<n*n; i++)
    //{
    //    printf("%d ", Ai[i]);
    //}
    //printf("\n");
    
    umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null) ;
    umfpack_di_free_symbolic (&Symbolic);
    umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null);
    umfpack_di_free_numeric (&Numeric) ;

    //for(i=0;i<n;i++) 
    //    printf("x[%d]=%g\n", i, x[i]) ;
    return 0;
}

 

后面会测试那种速度较快。