UMFPACK使用调用（三）

在调用UMFPACK的过程中，只需要关心Ap Ai Ax的产生，通过Eigen库，先让矩阵A以稀疏矩阵格式存储（知道矩阵A的非零元素的分布），调用UMFPACK成功

//#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=90000;

    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()]<<" ";
    double  *b=new double[n];
    for(int i=0;i<n;i++)
        b[i]=1;
    double *x=new double [n];

    double *null =(double *)NULL ;
    void *Symbolic, *Numeric ;

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

    return 0;
}

 参考：http://www.cnblogs.com/kmliang/archive/2013/03/14/2958852.html