portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Rework the linop support to get different forms of red black schur solver

Browse Source 
Moo on diag, or MooInv Moe MeeInv Meo
This commit is contained in:
Peter Boyle 2015-06-05 10:17:10 +01:00
parent 58a4f32298
commit 7f6304fac3
16 changed files with 155 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

154
lib/algorithms/LinearOperator.h
View File

@ -18,84 +18,144 @@ 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;
};
/////////////////////////////////////////////////////////////////////////////////////////////
// Hermitian operators are self adjoint and only require Op to be defined, so refine the base
/////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class HermitianOperatorBase : public LinearOperatorBase<Field> {
public:
virtual void OpAndNorm(const Field &in, Field &out,double &n1,double &n2)=0;
void AdjOp(const Field &in, Field &out) {
Op(in,out);
};
void Op(const Field &in, Field &out) {
double n1,n2;
OpAndNorm(in,out,n1,n2);
};
};
/////////////////////////////////////////////////////////////////////////////////////////////
// Whereas non hermitian takes a generic sparse matrix (e.g. lattice action)
// conforming to sparse matrix interface and builds the full checkerboard non-herm operator
// Op and AdjOp distinct.
// By sharing the class for Sparse Matrix across multiple operator wrappers, we can share code
// between RB and non-RB variants. Sparse matrix is like the fermion action def, and then
// the wrappers implement the specialisation of "Op" and "AdjOp" to the cases minimising
// replication of code.
//
// I'm not entirely happy with implementation; to share the Schur code between herm and non-herm
// while still having a "OpAndNorm" in the abstract base I had to implement it in both cases
// with an assert trap in the non-herm. This isn't right; there must be a better C++ way to
// do it, but I fear it required multiple inheritance and mixed in abstract base classes
/////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
// Construct herm op from non-herm matrix
////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class NonHermitianOperator : public LinearOperatorBase<Field> {
class MdagMLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
public:
NonHermitianOperator(Matrix &Mat): _Mat(Mat){};
MdagMLinearOperator(Matrix &Mat): _Mat(Mat){};
void Op (const Field &in, Field &out){
_Mat.M(in,out);
}
void AdjOp (const Field &in, Field &out){
_Mat.Mdag(in,out);
}
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
_Mat.MdagM(in,out,n1,n2);
}
};
////////////////////////////////////////////////////////////////////////////////////
// Redblack Non hermitian wrapper
////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
// Wrap an already herm matrix
////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class NonHermitianCheckerBoardedOperator : public LinearOperatorBase<Field> {
class HermitianLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
public:
NonHermitianCheckerBoardedOperator(Matrix &Mat): _Mat(Mat){};
HermitianLinearOperator(Matrix &Mat): _Mat(Mat){};
void Op (const Field &in, Field &out){
_Mat.Mpc(in,out);
_Mat.M(in,out);
}
void AdjOp (const Field &in, Field &out){ //
_Mat.MpcDag(in,out);
void AdjOp (const Field &in, Field &out){
_Mat.M(in,out);
}
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
ComplexD dot;
_Mat.M(in,out);
dot= innerProduct(in,out);
n1=real(dot);
dot = innerProduct(out,out);
n2=real(dot);
}
};
////////////////////////////////////////////////////////////////////////////////////
// Hermitian wrapper
////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class HermitianOperator : public HermitianOperatorBase<Field> {
Matrix &_Mat;
//////////////////////////////////////////////////////////
// Even Odd Schur decomp operators; there are several
// ways to introduce the even odd checkerboarding
//////////////////////////////////////////////////////////
template<class Field>
class SchurOperatorBase : public LinearOperatorBase<Field> {
public:
HermitianOperator(Matrix &Mat): _Mat(Mat) {};
void OpAndNorm(const Field &in, Field &out,double &n1,double &n2){
return _Mat.MdagM(in,out,n1,n2);
virtual RealD Mpc (const Field &in, Field &out) =0;
virtual RealD MpcDag (const Field &in, Field &out) =0;
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
Field tmp(in._grid);
ni=Mpc(in,tmp);
no=MpcDag(tmp,out);
}
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
MpcDagMpc(in,out,n1,n2);
}
void Op (const Field &in, Field &out){
Mpc(in,out);
}
void AdjOp (const Field &in, Field &out){
MpcDag(in,out);
}
};
////////////////////////////////////////////////////////////////////////////////////
// Hermitian CheckerBoarded wrapper
////////////////////////////////////////////////////////////////////////////////////
template<class Matrix,class Field>
class HermitianCheckerBoardedOperator : public HermitianOperatorBase<Field> {
class SchurDiagMooeeOperator : public SchurOperatorBase<Field> {
Matrix &_Mat;
public:
HermitianCheckerBoardedOperator(Matrix &Mat): _Mat(Mat) {};
void OpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
_Mat.MpcDagMpc(in,out,n1,n2);
SchurDiagMooeeOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in._grid);
_Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
_Mat.Mooee(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in._grid);
_Mat.MeooeDag(in,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeDag(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
};
template<class Matrix,class Field>
class SchurDiagOneOperator : public SchurOperatorBase<Field> {
Matrix &_Mat;
public:
SchurDiagOneOperator (Matrix &Mat): _Mat(Mat){};
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in._grid);
_Mat.Meooe(in,tmp);
_Mat.MooeeInv(tmp,out);
_Mat.Meooe(out,tmp);
_Mat.MooeeInv(tmp,out);
return axpy_norm(out,-1.0,tmp,in);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in._grid);
_Mat.MooeeInvDag(in,out);
_Mat.MeooeDag(out,tmp);
_Mat.MooeeInvDag(tmp,out);
_Mat.MeooeDag(out,tmp);
return axpy_norm(out,-1.0,tmp,in);
}
};
@ -106,10 +166,6 @@ namespace Grid {
public:
virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) = 0;
};
template<class Field> class HermitianOperatorFunction {
public:
virtual void operator() (HermitianOperatorBase<Field> &Linop, const Field &in, Field &out) = 0;
};
// FIXME : To think about

26
lib/algorithms/SparseMatrix.h
View File

@ -36,32 +36,6 @@ namespace Grid {
virtual void MooeeDag (const Field &in, Field &out)=0;
virtual void MooeeInvDag (const Field &in, Field &out)=0;
// Schur decomp operators
virtual RealD Mpc (const Field &in, Field &out) {
Field tmp(in._grid);
Meooe(in,tmp);
MooeeInv(tmp,out);
Meooe(out,tmp);
Mooee(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
virtual RealD MpcDag (const Field &in, Field &out){
Field tmp(in._grid);
MeooeDag(in,tmp);
MooeeInvDag(tmp,out);
MeooeDag(out,tmp);
MooeeDag(in,out);
return axpy_norm(out,-1.0,tmp,out);
}
virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) {
Field tmp(in._grid);
ni=Mpc(in,tmp);
no=MpcDag(tmp,out);
}
};
}

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

@ -9,17 +9,17 @@ namespace Grid {
/////////////////////////////////////////////////////////////
template<class Field>
class ConjugateGradient : public HermitianOperatorFunction<Field> {
class ConjugateGradient : public OperatorFunction<Field> {
public:
RealD Tolerance;
Integer MaxIterations;
int verbose;
ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) {
verbose=0;
verbose=1;
};
void operator() (HermitianOperatorBase<Field> &Linop,const Field &src, Field &psi){
void operator() (LinearOperatorBase<Field> &Linop,const Field &src, Field &psi){
psi.checkerboard = src.checkerboard;
conformable(psi,src);
@ -33,7 +33,7 @@ public:
//Initial residual computation & set up
RealD guess = norm2(psi);
Linop.OpAndNorm(psi,mmp,d,b);
Linop.HermOpAndNorm(psi,mmp,d,b);
r= src-mmp;
p= r;
@ -65,7 +65,7 @@ public:
c=cp;
Linop.OpAndNorm(p,mmp,d,qq);
Linop.HermOpAndNorm(p,mmp,d,qq);
RealD qqck = norm2(mmp);
ComplexD dck = innerProduct(p,mmp);
@ -86,19 +86,10 @@ public:
if (verbose) std::cout<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
// Hack
if (0) {
Field tt(src);
Linop.Op(psi,mmp);
tt=mmp-src;
RealD resnorm = norm2(tt);
std::cout<<"ConjugateGradient: Iteration " <<k<<" true residual "<<resnorm << " computed " << cp <<std::endl;
}
// Stopping condition
if ( cp <= rsq ) {
Linop.Op(psi,mmp);
Linop.HermOpAndNorm(psi,mmp,d,qq);
p=mmp-src;
RealD mmpnorm = sqrt(norm2(mmp));

15
lib/algorithms/iterative/SchurRedBlack.h
View File

@ -40,16 +40,16 @@ namespace Grid {
// Take a matrix and form a Red Black solver calling a Herm solver
// Use of RB info prevents making SchurRedBlackSolve conform to standard interface
///////////////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class SchurRedBlackSolve {
template<class Field> class SchurRedBlackDiagMooeeSolve {
private:
HermitianOperatorFunction<Field> & _HermitianRBSolver;
OperatorFunction<Field> & _HermitianRBSolver;
int CBfactorise;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick
/////////////////////////////////////////////////////
SchurRedBlackSolve(HermitianOperatorFunction<Field> &HermitianRBSolver) :
SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver) :
_HermitianRBSolver(HermitianRBSolver)
{
CBfactorise=0;
@ -62,6 +62,8 @@ namespace Grid {
// FIXME use CBfactorise to control schur decomp
GridBase *grid = _Matrix.RedBlackGrid();
GridBase *fgrid= _Matrix.Grid();
SchurDiagMooeeOperator<Matrix,Field> _HermOpEO(_Matrix);
Field src_e(grid);
Field src_o(grid);
@ -80,12 +82,13 @@ namespace Grid {
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd);
_Matrix.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd);
// get the right MpcDag
_HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd);
//////////////////////////////////////////////////////////////
// Call the red-black solver
//////////////////////////////////////////////////////////////
HermitianCheckerBoardedOperator<Matrix,Field> _HermOpEO(_Matrix);
std::cout << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
@ -105,7 +108,7 @@ namespace Grid {
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout << "SchurRedBlack solver true unprec resid "<< sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
std::cout << "SchurRedBlackDiagMooee solver true unprec resid "<< sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
}
};

6
tests/Test_cayley_cg.cc
View File

@ -130,7 +130,7 @@ void TestCGunprec(What & Ddwf,
LatticeFermion src (FGrid); random(*RNG5,src);
LatticeFermion result(FGrid); result=zero;
HermitianOperator<What,LatticeFermion> HermOp(Ddwf);
MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);
@ -149,7 +149,7 @@ void TestCGprec(What & Ddwf,
pickCheckerboard(Odd,src_o,src);
result_o=zero;
HermitianCheckerBoardedOperator<What,LatticeFermion> HermOpEO(Ddwf);
SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOpEO,src_o,result_o);
}
@ -167,6 +167,6 @@ void TestCGschur(What & Ddwf,
LatticeFermion result(FGrid); result=zero;
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
SchurSolver(Ddwf,src,result);
}

11
tests/Test_cayley_even_odd.cc
View File

@ -214,7 +214,7 @@ void TestWhat(What & Ddwf,
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"= Test DiagMoo MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
random(*RNG5,phi);
@ -225,11 +225,12 @@ void TestWhat(What & Ddwf,
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
Ddwf.MpcDagMpc(chi_e,dchi_e,t1,t2);
Ddwf.MpcDagMpc(chi_o,dchi_o,t1,t2);
SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
Ddwf.MpcDagMpc(phi_e,dphi_e,t1,t2);
Ddwf.MpcDagMpc(phi_o,dphi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);

6
tests/Test_contfrac_cg.cc
View File

@ -121,7 +121,7 @@ void TestCGunprec(What & Ddwf,
LatticeFermion src (FGrid); random(*RNG5,src);
LatticeFermion result(FGrid); result=zero;
HermitianOperator<What,LatticeFermion> HermOp(Ddwf);
MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);
@ -140,7 +140,7 @@ void TestCGprec(What & Ddwf,
pickCheckerboard(Odd,src_o,src);
result_o=zero;
HermitianCheckerBoardedOperator<What,LatticeFermion> HermOpEO(Ddwf);
SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOpEO,src_o,result_o);
}
@ -158,6 +158,6 @@ void TestCGschur(What & Ddwf,
LatticeFermion result(FGrid); result=zero;
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
SchurSolver(Ddwf,src,result);
}

11
tests/Test_contfrac_even_odd.cc
View File

@ -211,11 +211,12 @@ void TestWhat(What & Ddwf,
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
Ddwf.MpcDagMpc(chi_e,dchi_e,t1,t2);
Ddwf.MpcDagMpc(chi_o,dchi_o,t1,t2);
Ddwf.MpcDagMpc(phi_e,dphi_e,t1,t2);
Ddwf.MpcDagMpc(phi_o,dphi_o,t1,t2);
SchurDiagMooeeOperator<What,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);

2
tests/Test_dwf_cg_prec.cc
View File

@ -50,7 +50,7 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd,src_o,src);
result_o=zero;
HermitianCheckerBoardedOperator<DomainWallFermion,LatticeFermion> HermOpEO(Ddwf);
SchurDiagMooeeOperator<DomainWallFermion,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOpEO,src_o,result_o);

2
tests/Test_dwf_cg_schur.cc
View File

@ -46,7 +46,7 @@ int main (int argc, char ** argv)
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
SchurSolver(Ddwf,src,result);
Grid_finalize();

2
tests/Test_dwf_cg_unprec.cc
View File

@ -45,7 +45,7 @@ int main (int argc, char ** argv)
RealD M5=1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
HermitianOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);

10
tests/Test_dwf_even_odd.cc
View File

@ -186,11 +186,13 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
Ddwf.MpcDagMpc(chi_e,dchi_e,t1,t2);
Ddwf.MpcDagMpc(chi_o,dchi_o,t1,t2);
Ddwf.MpcDagMpc(phi_e,dphi_e,t1,t2);
Ddwf.MpcDagMpc(phi_o,dphi_o,t1,t2);
SchurDiagMooeeOperator<DomainWallFermion,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);

2
tests/Test_wilson_cg_prec.cc
View File

@ -53,7 +53,7 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd,src_o,src);
result_o=zero;
HermitianCheckerBoardedOperator<WilsonFermion,LatticeFermion> HermOpEO(Dw);
SchurDiagMooeeOperator<WilsonFermion,LatticeFermion> HermOpEO(Dw);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOpEO,src_o,result_o);

2
tests/Test_wilson_cg_schur.cc
View File

@ -40,7 +40,7 @@ int main (int argc, char ** argv)
WilsonFermion Dw(Umu,Grid,RBGrid,mass);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
SchurSolver(Dw,src,result);

2
tests/Test_wilson_cg_unprec.cc
View File

@ -49,7 +49,7 @@ int main (int argc, char ** argv)
RealD mass=0.5;
WilsonFermion Dw(Umu,Grid,RBGrid,mass);
HermitianOperator<WilsonFermion,LatticeFermion> HermOp(Dw);
MdagMLinearOperator<WilsonFermion,LatticeFermion> HermOp(Dw);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);

9
tests/Test_wilson_even_odd.cc
View File

@ -177,11 +177,12 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
Dw.MpcDagMpc(chi_e,dchi_e,t1,t2);
Dw.MpcDagMpc(chi_o,dchi_o,t1,t2);
SchurDiagMooeeOperator<WilsonFermion,LatticeFermion> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
Dw.MpcDagMpc(phi_e,dphi_e,t1,t2);
Dw.MpcDagMpc(phi_o,dphi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);