下面用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 
|
本类阅读TOP10