1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-14 09:45:36 +00:00
Grid/lib/algorithms/iterative/ConjugateGradient.h

106 lines
3.1 KiB
C
Raw Normal View History

2015-05-17 00:19:03 +01:00
#ifndef GRID_CONJUGATE_GRADIENT_H
#define GRID_CONJUGATE_GRADIENT_H
namespace Grid {
/////////////////////////////////////////////////////////////
// Base classes for iterative processes based on operators
// single input vec, single output vec.
/////////////////////////////////////////////////////////////
template<class Field>
class ConjugateGradient : public OperatorFunction<Field> {
public:
RealD Tolerance;
Integer MaxIterations;
ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
};
2015-05-17 00:19:03 +01:00
void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi){
2015-05-17 00:19:03 +01:00
psi.checkerboard = src.checkerboard;
conformable(psi,src);
2015-05-17 00:19:03 +01:00
RealD cp,c,a,d,b,ssq,qq,b_pred;
Field p(src);
Field mmp(src);
Field r(src);
//Initial residual computation & set up
RealD guess = norm2(psi);
Linop.HermOpAndNorm(psi,mmp,d,b);
2015-05-17 00:19:03 +01:00
r= src-mmp;
p= r;
a =norm2(p);
cp =a;
ssq=norm2(src);
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: cp,r "<<cp <<std::endl;
std::cout<<GridLogIterative <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
RealD rsq = Tolerance* Tolerance*ssq;
2015-05-17 00:19:03 +01:00
//Check if guess is really REALLY good :)
if ( cp <= rsq ) {
return;
2015-05-17 00:19:03 +01:00
}
std::cout<<GridLogIterative << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
2015-05-17 00:19:03 +01:00
int k;
for (k=1;k<=MaxIterations;k++){
2015-05-17 00:19:03 +01:00
c=cp;
Linop.HermOpAndNorm(p,mmp,d,qq);
RealD qqck = norm2(mmp);
ComplexD dck = innerProduct(p,mmp);
2015-05-17 00:19:03 +01:00
a = c/d;
b_pred = a*(a*qq-d)/c;
cp = axpy_norm(r,-a,mmp,r);
2015-05-17 00:19:03 +01:00
b = cp/c;
// Fuse these loops ; should be really easy
psi= a*p+psi;
p = p*b+r;
std::cout<<GridLogIterative<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
2015-05-17 00:19:03 +01:00
// Stopping condition
if ( cp <= rsq ) {
Linop.HermOpAndNorm(psi,mmp,d,qq);
2015-05-17 00:19:03 +01:00
p=mmp-src;
RealD mmpnorm = sqrt(norm2(mmp));
RealD psinorm = sqrt(norm2(psi));
RealD srcnorm = sqrt(norm2(src));
RealD resnorm = sqrt(norm2(p));
RealD true_residual = resnorm/srcnorm;
std::cout<<GridLogMessage<<"ConjugateGradient: Converged on iteration " <<k
<<" computed residual "<<sqrt(cp/ssq)
<<" true residual "<<true_residual
<<" target "<<Tolerance<<std::endl;
return;
2015-05-17 00:19:03 +01:00
}
}
std::cout<<GridLogMessage<<"ConjugateGradient did NOT converge"<<std::endl;
assert(0);
2015-05-17 00:19:03 +01:00
}
};
}
#endif