zoukankan      html  css  js  c++  java
  • UMFPACK使用调用(三)

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

    View Code
    //#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

  • 相关阅读:
    Populating Next Right Pointers in Each Node II
    Populating Next Right Pointers in Each Node
    Construct Binary Tree from Preorder and Inorder Traversal
    Construct Binary Tree from Inorder and Postorder Traversal
    Path Sum
    Symmetric Tree
    Solve Tree Problems Recursively
    632. Smallest Range(priority_queue)
    609. Find Duplicate File in System
    poj3159最短路spfa+邻接表
  • 原文地址:https://www.cnblogs.com/kmliang/p/2962553.html
Copyright © 2011-2022 走看看