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

Got unpreconditioned conjugate gradient to run and converge on a random (uniform random,

Browse Source 
not even SU(3) for now) gauge field. Convergence history is correctly indepdendent of decomposition
on 1,2,4,8,16 mpi tasks.
Found a couple of simd bugs which required fixed and enhanced the Grid_simd.cc test suite.
Implemented the Mdag, M, MdagM, Meooe Mooee schur type stuff in the wilson dop.
This commit is contained in:
Peter Boyle
2015-05-19 13:57:35 +01:00
parent 05d862782f
commit 4dba8522a1
33 changed files with 566 additions and 316 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
lib/algorithms/LinearOperator.h
View File

@ -7,8 +7,8 @@ namespace Grid {
// LinearOperators Take a something and return a something.
/////////////////////////////////////////////////////////////////////////////////////////////
//
// Hopefully linearity is satisfied and the AdjOp is indeed the Hermitian conjugate (transpose if real):
//
// Hopefully linearity is satisfied and the AdjOp is indeed the Hermitian conjugateugate (transpose if real):
//SBase
// i) F(a x + b y) = aF(x) + b F(y).
// ii) <x|Op|y> = <y|AdjOp|x>^\ast
//
@ -25,12 +25,13 @@ namespace Grid {
/////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class HermitianOperatorBase : public LinearOperatorBase<Field> {
public:
virtual RealD OpAndNorm(const Field &in, Field &out);
virtual void OpAndNorm(const Field &in, Field &out,double &n1,double &n2);
void AdjOp(const Field &in, Field &out) {
Op(in,out);
};
void Op(const Field &in, Field &out) {
OpAndNorm(in,out);
double n1,n2;
OpAndNorm(in,out,n1,n2);
};
};
@ -80,8 +81,8 @@ namespace Grid {
Matrix &_Mat;
public:
HermitianOperator(Matrix &Mat): _Mat(Mat) {};
RealD OpAndNorm(const Field &in, Field &out){
return _Mat.MdagM(in,out);
void OpAndNorm(const Field &in, Field &out,double &n1,double &n2){
return _Mat.MdagM(in,out,n1,n2);
}
};

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

@ -18,6 +18,7 @@ public:
std::cout << Tolerance<<std::endl;
};
void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi) {assert(0);};
void operator() (HermitianOperatorBase<Field> &Linop,const Field &src, Field &psi){
RealD cp,c,a,d,b,ssq,qq,b_pred;
@ -61,21 +62,27 @@ public:
Linop.OpAndNorm(p,mmp,d,qq);
// std::cout <<std::setprecision(4)<< "ConjugateGradient: d,qq "<<d<< " "<<qq <<std::endl;
a = c/d;
b_pred = a*(a*qq-d)/c;
cp = axpy_norm(r,mmp,r,-a);
// std::cout <<std::setprecision(4)<< "ConjugateGradient: a,bp "<<a<< " "<<b_pred <<std::endl;
cp = axpy_norm(r,-a,mmp,r);
b = cp/c;
// std::cout <<std::setprecision(4)<< "ConjugateGradient: cp,b "<<cp<< " "<<b <<std::endl;
// Fuse these loops ; should be really easy
psi= a*p+psi;
p = p*b+r;
std::cout << "Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
std::cout<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
// Stopping condition
if ( cp <= rsq ) {
Linop.Op(p,mmp);
Linop.Op(psi,mmp);
p=mmp-src;
RealD mmpnorm = sqrt(norm2(mmp));
@ -83,8 +90,11 @@ public:
RealD srcnorm = sqrt(norm2(src));
RealD resnorm = sqrt(norm2(p));
RealD true_residual = resnorm/srcnorm;
std::cout<<"ConjugateGradient: Converged on iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
std::cout<<"ConjugateGradient: true residual is "<<true_residual<<" sol "<<psinorm<<" src "<<srcnorm<<std::endl;
std::cout<<"ConjugateGradient: target residual was "<<Tolerance<<std::endl;
return;
}
}
std::cout<<"ConjugateGradient did NOT converge"<<std::endl;