portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-13 12:47:05 +01:00
Code Releases Activity

Getting closer to having a wilson solver... introducing a first and untested

Browse Source 
cut at Conjugate gradient. Also copied in Remez, Zolotarev, Chebyshev from
Mike Clark, Tony Kennedy and my BFM package respectively since we know we will
need these. I wanted the structure of

algorithms/approx
algorithms/iterative

etc.. to start taking shape.
This commit is contained in:
Peter Boyle
2015-05-18 07:47:05 +01:00
parent 7992346190
commit 11cb3e9a01
22 changed files with 1798 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files

49
lib/algorithms/iterative/ConjugateGradient.h
View File

@ -1,19 +1,25 @@
#ifndef GRID_CONJUGATE_GRADIENT_H
#define GRID_CONJUGATE_GRADIENT_H
#include<Grid.h>
#include<algorithms/LinearOperator.h>
namespace Grid {
template<class Field> class ConjugateGradient : public IterativeProcess<Field> {
public:
/////////////////////////////////////////////////////////////
// Base classes for iterative processes based on operators
// single input vec, single output vec.
/////////////////////////////////////////////////////////////
ConjugateGradient(RealD tol,Integer maxit): IterativeProces(tol,maxit) {};
template<class Field>
class ConjugateGradient : public OperatorFunction<Field> {
public:
RealD Tolerance;
Integer MaxIterations;
void operator() (HermitianOperatorBase<Field> *Linop,const Field &src, Field &psi){
ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
std::cout << Tolerance<<std::endl;
};
void operator() (HermitianOperatorBase<Field> &Linop,const Field &src, Field &psi){
RealD residual = Tolerance;
RealD cp,c,a,d,b,ssq,qq,b_pred;
Field p(src);
@ -32,24 +38,24 @@ namespace Grid {
cp =a;
ssq=norm2(src);
std::cout <<setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
std::cout <<setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
std::cout <<setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
std::cout <<setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
std::cout <<setprecision(4)<< "ConjugateGradient: r "<<cp <<std::endl;
std::cout <<setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: r "<<cp <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
RealD rsq = residual* residual*ssq;
RealD rsq = Tolerance* Tolerance*ssq;
//Check if guess is really REALLY good :)
if ( cp <= rsq ) {
return 0;
return;
}
std::cout << setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
std::cout << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
int k;
for (k=1;k<=max_iter;k++){
for (k=1;k<=MaxIterations;k++){
c=cp;
@ -62,10 +68,10 @@ namespace Grid {
b = cp/c;
// Fuse these loops ; should be really easy
// Update psi; New (conjugate/M-orthogonal) search direction
psi= a*p+psi;
p = p*b+r;
std::cout << "Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
// Stopping condition
if ( cp <= rsq ) {
@ -78,10 +84,11 @@ namespace Grid {
RealD resnorm = sqrt(norm2(p));
RealD true_residual = resnorm/srcnorm;
std::cout<<"ConjugateGradient: true residual is "<<true_residual<<" sol "<<psinorm<<" src "<<srcnorm<<std::endl;
std::cout<<"ConjugateGradient: target residual was "<<residual<<std::endl;
return k;
std::cout<<"ConjugateGradient: target residual was "<<Tolerance<<std::endl;
}
}
std::cout<<"ConjugateGradient did NOT converge"<<std::endl;
assert(0);
}
};
}