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

Conjugate residual added

Browse Source
This commit is contained in:
Peter Boyle
2015-06-05 18:16:25 +01:00
parent 351c2905f5
commit 1a05882d7c
15 changed files with 181 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
lib/algorithms/LinearOperator.h
View File

@ -18,7 +18,8 @@ namespace Grid {
public:
virtual void Op (const Field &in, Field &out) = 0; // Abstract base
virtual void AdjOp (const Field &in, Field &out) = 0; // Abstract base
virtual void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2)=0;
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2)=0;
virtual void HermOp(const Field &in, Field &out)=0;
};
@ -48,9 +49,13 @@ namespace Grid {
void AdjOp (const Field &in, Field &out){
_Mat.Mdag(in,out);
}
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.MdagM(in,out,n1,n2);
}
void HermOp(const Field &in, Field &out){
RealD n1,n2;
HermOpAndNorm(in,out,n1,n2);
}
};
////////////////////////////////////////////////////////////////////
@ -67,7 +72,7 @@ namespace Grid {
void AdjOp (const Field &in, Field &out){
_Mat.M(in,out);
}
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
ComplexD dot;
_Mat.M(in,out);
@ -78,6 +83,9 @@ namespace Grid {
dot = innerProduct(out,out);
n2=real(dot);
}
void HermOp(const Field &in, Field &out){
_Mat.M(in,out);
}
};
//////////////////////////////////////////////////////////
@ -98,6 +106,10 @@ namespace Grid {
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
MpcDagMpc(in,out,n1,n2);
}
void HermOp(const Field &in, Field &out){
RealD n1,n2;
HermOpAndNorm(in,out,n1,n2);
}
void Op (const Field &in, Field &out){
Mpc(in,out);
}

1
lib/algorithms/SparseMatrix.h
View File

@ -18,6 +18,7 @@ namespace Grid {
Field tmp (in._grid);
ni=M(in,tmp);
no=Mdag(tmp,out);
std::cout << "MdagM "<< ni<<" "<<no<<std::endl;
}
};

85
lib/algorithms/iterative/ConjugateResidual.h Normal file
View File

@ -0,0 +1,85 @@
#ifndef GRID_CONJUGATE_RESIDUAL_H
#define GRID_CONJUGATE_RESIDUAL_H
namespace Grid {
/////////////////////////////////////////////////////////////
// Base classes for iterative processes based on operators
// single input vec, single output vec.
/////////////////////////////////////////////////////////////
template<class Field>
class ConjugateResidual : public OperatorFunction<Field> {
public:
RealD Tolerance;
Integer MaxIterations;
int verbose;
ConjugateResidual(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
verbose=1;
};
void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi){
RealD a, b, c, d;
RealD cp, ssq,rsq;
RealD rAr, rAAr, rArp;
RealD pAp, pAAp;
GridBase *grid = src._grid;
psi=zero;
Field r(grid), p(grid), Ap(grid), Ar(grid);
r=src;
p=src;
Linop.HermOpAndNorm(p,Ap,pAp,pAAp);
Linop.HermOpAndNorm(r,Ar,rAr,rAAr);
std::cout << "pAp, pAAp"<< pAp<<" "<<pAAp<<std::endl;
std::cout << "rAr, rAAr"<< rAr<<" "<<rAAr<<std::endl;
cp =norm2(r);
ssq=norm2(src);
rsq=Tolerance*Tolerance*ssq;
std::cout<<"ConjugateResidual: iteration " <<0<<" residual "<<cp<< " target"<< rsq<<std::endl;
for(int k=1;k<MaxIterations;k++){
a = rAr/pAAp;
axpy(psi,a,p,psi);
cp = axpy_norm(r,-a,Ap,r);
rArp=rAr;
Linop.HermOpAndNorm(r,Ar,rAr,rAAr);
b =rAr/rArp;
axpy(p,b,p,r);
pAAp=axpy_norm(Ap,b,Ap,Ar);
std::cout<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
if(cp<rsq) {
Linop.HermOp(psi,Ap);
axpy(r,-1.0,src,Ap);
RealD true_resid = norm2(r);
std::cout<<"ConjugateResidual: Converged on iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
std::cout<<"ConjugateResidual: true residual is "<<true_resid<<std::endl;
std::cout<<"ConjugateResidual: target residual was "<<Tolerance <<std::endl;
return;
}
}
std::cout<<"ConjugateResidual did NOT converge"<<std::endl;
assert(0);
}
};
}
#endif