下面用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) int maxit = 150; // maximum iteration CompCol_Mat_double Jacobi(nUnknown, nUnknown, nNonZero, val, row_ind, col_ptr); MATRIX_double H(restart+1, restart, 0.0); // storage for upper Hessenberg H; for(int i=0; i<nUnknown; i++) rhs(i) =i+1; /********************************************************************** cout << "GMRES flag = " << results << endl; |
运行结果:
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
用十万阶sparse matrix测试性能满意。