LA_library/conjgrad.h

97 lines
2.8 KiB
C
Raw Normal View History

2008-02-26 14:55:23 +01:00
/*
LA: linear algebra C++ interface library
Copyright (C) 2008 Jiri Pittner <jiri.pittner@jh-inst.cas.cz> or <jiri@pittnerovi.com>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
2005-02-04 15:31:42 +01:00
#include "vec.h"
#include "smat.h"
#include "mat.h"
#include "sparsemat.h"
#include "nonclass.h"
2006-10-23 20:20:58 +02:00
#include <iomanip>
2005-02-04 15:31:42 +01:00
//conjugate gradient solution of a linear system
2006-10-25 10:29:57 +02:00
//matrix can be any class which has nrows(), ncols(), diagonalof() and gemv() available
2005-02-04 15:31:42 +01:00
//does not even have to be explicitly stored
2006-10-23 20:20:58 +02:00
//Conjugate gradient algorithm, cf. Bulirsch-Stoer book
2005-02-04 15:31:42 +01:00
template<typename T, typename Matrix>
2008-12-18 17:07:06 +01:00
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)
2006-10-23 20:20:58 +02:00
{
int m=bigmat.nrows();
int n=bigmat.ncols();
if(x.size()!=n || b.size() != m) laerror("incompatible vectors and matrix sizes in conjgrad");
if(m!=n) issquare=0;
double t,tt,bscal,ascal;
NRVec<T> p,rr, *r;
NRVec<T> q(m),s(m);
if(issquare) r=&s; else r = new NRVec<T>(m);
if(doguess)
{
2006-10-25 10:29:57 +02:00
bigmat.gemv(0,x,'t',-1.,b); //x.gemv(0,bigmat,'t',-1.,b);
2006-10-23 20:20:58 +02:00
if(precondition) bigmat.diagonalof(x,true);
x.normalize();
}
2006-10-25 10:29:57 +02:00
bigmat.gemv(0,s,'n',-1.,x); //s.gemv(0,bigmat,'n',-1.,x);
2006-10-23 20:20:58 +02:00
s+=b;
2006-10-25 10:29:57 +02:00
if(!issquare) bigmat.gemv(0,*r,'t',1,s); //(*r).gemv(0,bigmat,'t',1,s);
2006-10-23 20:20:58 +02:00
rr= *r;
if(precondition) bigmat.diagonalof(rr,true);
p=rr;
for(int iter=0; iter<= itmax; iter++)
{
double err=p.norm();
2008-12-18 16:48:25 +01:00
if(verbose)
{
cout << "conjgrad: iter= "<<iter<<" err= "<<
2006-10-23 20:20:58 +02:00
setiosflags(ios::scientific)<<setprecision(8) <<err<<
resetiosflags(ios::scientific)<<setprecision(12)<<"\n";
2008-12-18 16:48:25 +01:00
cout.flush();
}
2006-10-27 15:40:09 +02:00
if(err <= tol)
{
if(!issquare) delete r;
return true;
}
2006-10-23 20:20:58 +02:00
2006-10-25 10:29:57 +02:00
bigmat.gemv(0,q,'n',1,p); //q.gemv(0,bigmat,'n',1,p);
2006-10-23 20:20:58 +02:00
tt= (*r) * rr;
t=issquare?p*q:q*q;
if(!t) {if(!issquare) delete r; laerror("conjgrad: singular matrix 1");}
ascal=tt/t;
x.axpy(ascal,p);
s.axpy(-ascal,q);
2006-10-25 10:29:57 +02:00
if(!issquare) bigmat.gemv(0,*r,'t',1,s); //(*r).gemv(0,bigmat,'t',1,s);
2006-10-23 20:20:58 +02:00
rr= *r;
if(precondition) bigmat.diagonalof(rr,true);
if(!tt) {if(!issquare) delete r; laerror("conjgrad: singular matrix 2");}
bscal= ((*r)*rr)/tt;
rr.axpy(bscal,p);
p=rr;
}
if(!issquare) delete r;
2006-10-27 15:40:09 +02:00
return false;
2006-10-23 20:20:58 +02:00
}