改写UMFPACK算例中的压缩方式

在UMFPACK的官方文档 UMFPACK Version 5.2.0 User Guide中的5.3节中的例子

稀疏矩阵A和右端向量b 为

 A=[2 3 0 0 0;3 0 4 0 6 ;0 -1 -3 2 0;0 0 1 0 0;0 4 2 0 1]
 b=[8.,45.,-3.,3.,19.]

其解为x=[1 2 3 4 5]'

这个算例的源代码为

View Code

 1 #include  <stdio.h>
 2 #include  "umfpack.h"
 3 int       n  =  5  ;
 4 int       Ap  [  ]  =  {0,  2,  5,  9,  10,  12}  ;
 5 int       Ai  [  ]  =  {  0,    1,    0,     2,    4,    1,    2,    3,     4,    2,    1,    4}  ;
 6 14double  Ax  [  ]  =  {2.,  3.,  3.,  -1.,  4.,  4.,  -3.,  1.,  2.,  2.,  6.,  1.}  ;
 7 double  b  [  ]  =  {8.,  45.,  -3.,  3.,  19.}  ;
 8 double  x  [5]  ;
 9 int  main  (void)
10 {
11     double  *null  =  (double  *)  NULL  ;
12     int  i  ;
13     void  *Symbolic,  *Numeric  ;
14     (void)  umfpack_di_symbolic  (n,  n,  Ap,  Ai,  Ax,  &Symbolic,  null,  null)  ;
15     (void)  umfpack_di_numeric  (Ap,  Ai,  Ax,  Symbolic,  &Numeric,  null,  null)  ;
16     umfpack_di_free_symbolic  (&Symbolic)  ;
17     (void)  umfpack_di_solve  (UMFPACK_A,  Ap,  Ai,  Ax,  x,  b,  Numeric,  null,  null)  ;
18     umfpack_di_free_numeric  (&Numeric)  ;
19     for  (i  =  0  ;  i  <  n  ;  i++)  printf  ("x  [%d]  =  %g\n",  i,  x  [i])  ;
20     return  (0)  ;
21 }

这个地方Ap, Ai, and Ax 需要手动输入，现改为自动输入

View Code

/* -------------------------------------------------------------------------- */
/* UMFPACK Copyright (c) Timothy A. Davis, CISE,                              */
/* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
/* -------------------------------------------------------------------------- */

#include <stdio.h>
#include "umfpack.h"
#include<vector>
#include<iostream>
using namespace std;

int main (void)
{
    int A[5][5]={{2,3,0,0,0},{3,0,4,0,6},{0,-1,-3,2,0},{0,0,1,0,0},{0,4,2,0,1}};
    int    n = 5 ;
    int    Ap [ ] = {0, 2, 5, 9, 10, 12} ;
    //int *Ap=new int [n+1];
    //Ap [ ]= {0, 2, 5, 9, 10, 12} ;
    int Apm[6]={0};//相当于Ap
    double epsilon=0.00001;
    int NZnum=0;//矩阵非零元的个数
    vector<int> Axm;//相当于Ax
    vector<int> Aim;//相当于Ai
    for(int j=0;j<5;j++)
    {
        for(int i=0;i<5;i++)
        {
            if(abs(A[i][j])>epsilon)
            {
                //Axm[NZnum]=A[i][j];printf("%d\n",Axm[NZnum]);
                Axm.push_back(A[i][j]);
                Aim.push_back(i);
                NZnum++;
            }
        }
        Apm[j+1]=NZnum;
        //printf("%d\n",Apm[j+1]);
    }
    for (int i=0; i<Axm.size(); i++)
    {
        printf("%d ", Axm[i]);
    }
    cout<<endl;
    for (int i=0; i<Aim.size(); i++)
    {
        printf("%d ", Aim[i]);
    }
    cout<<endl;
    int    Ai [ ] = { 0,  1,  0,   2,  4,  1,  2,  3,   4,  2,  1,  4} ;
    double Ax [ ] = {2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1.} ;
    double b [ ] = {8., 45., -3., 3., 19.} ;
    double x [5] ;
    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) ;
    for (i = 0 ; i < n ; i++) printf ("x [%d] = %g\n", i, x [i]) ;
    return (0) ;
}

可以达到同样目的采用的主要手段是vector