00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef _CONJGRAD_H_
00019 #define _CONJGRAD_H_
00020 #include "vec.h"
00021 #include "smat.h"
00022 #include "mat.h"
00023 #include "sparsemat.h"
00024 #include "nonclass.h"
00025 #include <iomanip>
00026
00027 namespace LA {
00028
00029
00030
00031
00032
00033
00034
00035
00036 template<typename T, typename Matrix>
00037 extern bool conjgrad(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=true, const double tol=1e-8, const int itmax=1000, const bool verbose=true, bool issquare=1,const bool precondition=1)
00038 {
00039 int m=bigmat.nrows();
00040 int n=bigmat.ncols();
00041
00042 if(x.size()!=n || b.size() != m) laerror("incompatible vectors and matrix sizes in conjgrad");
00043 if(m!=n) issquare=0;
00044
00045 double t,tt,bscal,ascal;
00046
00047 NRVec<T> p,rr, *r;
00048 NRVec<T> q(m),s(m);
00049 if(issquare) r=&s; else r = new NRVec<T>(m);
00050
00051 if(doguess)
00052 {
00053 bigmat.gemv(0,x,'t',-1.,b);
00054 if(precondition) bigmat.diagonalof(x,true);
00055 x.normalize();
00056 }
00057
00058 bigmat.gemv(0,s,'n',-1.,x);
00059 s+=b;
00060 if(!issquare) bigmat.gemv(0,*r,'t',1,s);
00061 rr= *r;
00062 if(precondition) bigmat.diagonalof(rr,true);
00063 p=rr;
00064
00065 for(int iter=0; iter<= itmax; iter++)
00066 {
00067 double err=p.norm();
00068 if(verbose)
00069 {
00070 std::cout << "conjgrad: iter= "<<iter<<" err= "<<
00071 std::setiosflags(std::ios::scientific)<<std::setprecision(8) <<err<<
00072 std::resetiosflags(std::ios::scientific)<<std::setprecision(12)<<"\n";
00073 std::cout.flush();
00074 }
00075 if(err <= tol)
00076 {
00077 if(!issquare) delete r;
00078 return true;
00079 }
00080
00081 bigmat.gemv(0,q,'n',1,p);
00082 tt= (*r) * rr;
00083 t=issquare?p*q:q*q;
00084 if(!t) {if(!issquare) delete r; laerror("conjgrad: singular matrix 1");}
00085 ascal=tt/t;
00086 x.axpy(ascal,p);
00087 s.axpy(-ascal,q);
00088 if(!issquare) bigmat.gemv(0,*r,'t',1,s);
00089 rr= *r;
00090 if(precondition) bigmat.diagonalof(rr,true);
00091 if(!tt) {if(!issquare) delete r; laerror("conjgrad: singular matrix 2");}
00092 bscal= ((*r)*rr)/tt;
00093 rr.axpy(bscal,p);
00094 p=rr;
00095 }
00096
00097 if(!issquare) delete r;
00098 return false;
00099 }
00100
00101 }
00102 #endif