用sparselib库解稀疏矩阵线性方程组

下面用GMRES（Generalized Minimum Residual Method）
演示用sparselib解线性方程组。
在matlab里可以用以下的命令，
GMRES(A,B,RESTART,TOL,MAXIT)；
类似得在vc++中可以这样：
#include <cstdlib>
#include <iostream>
#include "compcol_double.h"
#include "mvvtp.h"
#include "mvblasd.h"
#include "ilupre_double.h"
#include "gmres.h"
#include "spblas.h"
#include "mvm.h"
//#include MATRIX_H
//using namespace std;

int main(void)
{
   double val[] = {10, 3, 3, 9, 7, 8, 4, 9, 8, 7, 7, 9, -2, 5, 9, 2, 3, 13, -1};
   int row_ind[] = {0, 1, 3, 1, 2, 4, 5, 2, 3, 2, 3, 4,  0, 3, 4, 5, 1, 4,  5};
   int col_ptr[] = {0, 3, 7, 9, 12, 16, 19};

   int maxit = 150;                     // maximum iteration
   int nUnknown = 6;                    // unknown, the size of Jacobi
   int nNonZero = 19;                   // nonZero values in the matrix
   int results;
   int restart = 10;                    // restart iterations
   double tol = 1.e-6;                  // convergence tolerance

   CompCol_Mat_double  Jacobi(nUnknown, nUnknown, nNonZero, val, row_ind, col_ptr);
   //cout << Jacobi;
   CompCol_ILUPreconditioner_double M(Jacobi); // construct preconditioner

   MATRIX_double H(restart+1, restart, 0.0);   // storage for upper Hessenberg H;
   VECTOR_double xi(nUnknown, 0);
   VECTOR_double rhs(nUnknown);

   for(int i=0; i<nUnknown; i++)     rhs(i) =i+1;

/**********************************************************************
*      maxit AND tol WILL BE CHANGED AFTER ONE CALL OF GMRES,        **
*      SO FOR NEXT CALL, YOU SHOULD RESTORE THE OLD VALUE OF THEM    **
***********************************************************************/
   results = GMRES(Jacobi, xi, rhs, M, H, restart, maxit, tol); // call solver

   cout << "GMRES flag = " << results << endl;
   cout << "Iterations performed: " << maxit << endl;
   cout << "Tolerance achieved :" << tol << endl;
   for (i = 0; i < nUnknown; i ++){
     cout <<"xi["<<i<<"]="<<xi[i]<<"\n";
   }
   return results;
}

运行结果：
xi[0]=  0.248096
xi[1]=  0.705373
xi[2]=-1.49092
xi[3]= 1.64009
xi[4]= 0.740481
xi[5]=-1.69755
注意这里的稀疏矩阵用Harwell-Boeing格式存储,
    10     0     0     0    -2     0
     3     9     0     0     0     3
     0     7     9     7     0     0
     3     0     8     7     5     0
     0     8     0     9     9    13
     0     4     0     0     2    -1

具体可以参考
I.S. Duff,R.G.Grimes,and J.G.Lewis,Sparse matrix test problems,ACM Trans.Math.Soft

本文地址：http://com.8s8s.com/it/it23980.htm