mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-15 02:05:37 +00:00
Re-added Stag preconditioned opertor and preconditioned solvers, as the current versions are not consistent with CPS or MILC conventions.
This commit is contained in:
Chulwoo Jung
parent
5d44346be3
commit
0cccb3dc82
@ -8,6 +8,7 @@
|
|||||||
|
|
||||||
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||||
|
Author: Chulwoo Jung <chulwoo@bnl.gov>
|
||||||
|
|
||||||
This program is free software; you can redistribute it and/or modify
|
This program is free software; you can redistribute it and/or modify
|
||||||
it under the terms of the GNU General Public License as published by
|
it under the terms of the GNU General Public License as published by
|
||||||
@ -330,7 +331,61 @@ namespace Grid {
|
|||||||
assert(0);// Never need with staggered
|
assert(0);// Never need with staggered
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
|
// template<class Matrix,class Field> using SchurStagOperator = SchurStaggeredOperator<Matrix,Field>;
|
||||||
|
template<class Matrix,class Field>
|
||||||
|
// class SchurStagOperator : public LinearOperatorBase<Field> {
|
||||||
|
class SchurStagOperator : public SchurOperatorBase<Field> {
|
||||||
|
protected:
|
||||||
|
Matrix &_Mat;
|
||||||
|
public:
|
||||||
|
SchurStagOperator (Matrix &Mat): _Mat(Mat){};
|
||||||
|
virtual RealD Mpc (const Field &in, Field &out) {
|
||||||
|
Field tmp(in._grid);
|
||||||
|
Field tmp2(in._grid);
|
||||||
|
|
||||||
|
_Mat.Mooee(in,out);
|
||||||
|
_Mat.Mooee(out,tmp);
|
||||||
|
|
||||||
|
_Mat.Meooe(in,out);
|
||||||
|
_Mat.Meooe(out,tmp2);
|
||||||
|
|
||||||
|
return axpy_norm(out,-1.0,tmp2,tmp);
|
||||||
|
}
|
||||||
|
virtual RealD MpcDag (const Field &in, Field &out){
|
||||||
|
|
||||||
|
return Mpc(in,out);
|
||||||
|
}
|
||||||
|
#if 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);
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
|
||||||
|
n2 = Mpc(in,out);
|
||||||
|
ComplexD dot = innerProduct(in,out);
|
||||||
|
n1 = real(dot);
|
||||||
|
}
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
void AdjOp (const Field &in, Field &out){
|
||||||
|
MpcDag(in,out);
|
||||||
|
}
|
||||||
|
// Support for coarsening to a multigrid
|
||||||
|
void OpDiag (const Field &in, Field &out) {
|
||||||
|
assert(0); // must coarsen the unpreconditioned system
|
||||||
|
}
|
||||||
|
void OpDir (const Field &in, Field &out,int dir,int disp) {
|
||||||
|
assert(0);
|
||||||
|
}
|
||||||
|
|
||||||
|
};
|
||||||
|
|
||||||
#if 0
|
#if 0
|
||||||
// This is specific to (Z)mobius fermions
|
// This is specific to (Z)mobius fermions
|
||||||
|
|||||||
@ -164,7 +164,83 @@ namespace Grid {
|
|||||||
std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
|
std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
|
// template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
|
||||||
|
template<class Field> class SchurRedBlackStagSolve {
|
||||||
|
private:
|
||||||
|
OperatorFunction<Field> & _HermitianRBSolver;
|
||||||
|
int CBfactorise;
|
||||||
|
public:
|
||||||
|
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
// Wrap the usual normal equations Schur trick
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
SchurRedBlackStagSolve(OperatorFunction<Field> &HermitianRBSolver) :
|
||||||
|
_HermitianRBSolver(HermitianRBSolver)
|
||||||
|
{
|
||||||
|
CBfactorise=0;
|
||||||
|
};
|
||||||
|
|
||||||
|
template<class Matrix>
|
||||||
|
void operator() (Matrix & _Matrix,const Field &in, Field &out){
|
||||||
|
|
||||||
|
// FIXME CGdiagonalMee not implemented virtual function
|
||||||
|
// FIXME use CBfactorise to control schur decomp
|
||||||
|
GridBase *grid = _Matrix.RedBlackGrid();
|
||||||
|
GridBase *fgrid= _Matrix.Grid();
|
||||||
|
|
||||||
|
SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix);
|
||||||
|
|
||||||
|
Field src_e(grid);
|
||||||
|
Field src_o(grid);
|
||||||
|
Field sol_e(grid);
|
||||||
|
Field sol_o(grid);
|
||||||
|
Field tmp(grid);
|
||||||
|
Field Mtmp(grid);
|
||||||
|
Field resid(fgrid);
|
||||||
|
|
||||||
|
pickCheckerboard(Even,src_e,in);
|
||||||
|
pickCheckerboard(Odd ,src_o,in);
|
||||||
|
pickCheckerboard(Even,sol_e,out);
|
||||||
|
pickCheckerboard(Odd ,sol_o,out);
|
||||||
|
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
// src_o = Mdag * (source_o - Moe MeeInv source_e)
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
_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);
|
||||||
|
|
||||||
|
#if 0
|
||||||
|
// get the right MpcDag
|
||||||
|
// _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd);
|
||||||
|
#else
|
||||||
|
_Matrix.Mooee(tmp,src_o); assert(src_o.checkerboard ==Odd);
|
||||||
|
#endif
|
||||||
|
//////////////////////////////////////////////////////////////
|
||||||
|
// Call the red-black solver
|
||||||
|
//////////////////////////////////////////////////////////////
|
||||||
|
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
|
||||||
|
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////
|
||||||
|
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
|
||||||
|
///////////////////////////////////////////////////
|
||||||
|
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even);
|
||||||
|
src_e = src_e-tmp; assert( src_e.checkerboard ==Even);
|
||||||
|
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even);
|
||||||
|
|
||||||
|
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even);
|
||||||
|
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd );
|
||||||
|
|
||||||
|
// Verify the unprec residual
|
||||||
|
_Matrix.M(out,resid);
|
||||||
|
resid = resid-in;
|
||||||
|
RealD ns = norm2(in);
|
||||||
|
RealD nr = norm2(resid);
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage << "SchurRedBlackStag solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
///////////////////////////////////////////////////////////////////////////////////////////////////////
|
///////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||||
// Take a matrix and form a Red Black solver calling a Herm solver
|
// Take a matrix and form a Red Black solver calling a Herm solver
|
||||||
@ -403,5 +479,78 @@ namespace Grid {
|
|||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
template<class Field> class SchurRedBlackStagMixed {
|
||||||
|
private:
|
||||||
|
LinearFunction<Field> & _HermitianRBSolver;
|
||||||
|
int CBfactorise;
|
||||||
|
public:
|
||||||
|
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
// Wrap the usual normal equations Schur trick
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
SchurRedBlackStagMixed(LinearFunction<Field> &HermitianRBSolver) :
|
||||||
|
_HermitianRBSolver(HermitianRBSolver)
|
||||||
|
{
|
||||||
|
CBfactorise=0;
|
||||||
|
};
|
||||||
|
|
||||||
|
template<class Matrix>
|
||||||
|
void operator() (Matrix & _Matrix,const Field &in, Field &out){
|
||||||
|
|
||||||
|
// FIXME CGdiagonalMee not implemented virtual function
|
||||||
|
// FIXME use CBfactorise to control schur decomp
|
||||||
|
GridBase *grid = _Matrix.RedBlackGrid();
|
||||||
|
GridBase *fgrid= _Matrix.Grid();
|
||||||
|
|
||||||
|
SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix);
|
||||||
|
|
||||||
|
Field src_e(grid);
|
||||||
|
Field src_o(grid);
|
||||||
|
Field sol_e(grid);
|
||||||
|
Field sol_o(grid);
|
||||||
|
Field tmp(grid);
|
||||||
|
Field Mtmp(grid);
|
||||||
|
Field resid(fgrid);
|
||||||
|
|
||||||
|
pickCheckerboard(Even,src_e,in);
|
||||||
|
pickCheckerboard(Odd ,src_o,in);
|
||||||
|
pickCheckerboard(Even,sol_e,out);
|
||||||
|
pickCheckerboard(Odd ,sol_o,out);
|
||||||
|
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
// src_o = Mdag * (source_o - Moe MeeInv source_e)
|
||||||
|
/////////////////////////////////////////////////////
|
||||||
|
_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.Mooee(tmp,src_o); assert(src_o.checkerboard ==Odd);
|
||||||
|
//////////////////////////////////////////////////////////////
|
||||||
|
// Call the red-black solver
|
||||||
|
//////////////////////////////////////////////////////////////
|
||||||
|
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
|
||||||
|
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
|
||||||
|
_HermitianRBSolver(src_o,sol_o); assert(sol_o.checkerboard==Odd);
|
||||||
|
|
||||||
|
///////////////////////////////////////////////////
|
||||||
|
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
|
||||||
|
///////////////////////////////////////////////////
|
||||||
|
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even);
|
||||||
|
src_e = src_e-tmp; assert( src_e.checkerboard ==Even);
|
||||||
|
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even);
|
||||||
|
|
||||||
|
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even);
|
||||||
|
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd );
|
||||||
|
|
||||||
|
// Verify the unprec residual
|
||||||
|
_Matrix.M(out,resid);
|
||||||
|
resid = resid-in;
|
||||||
|
RealD ns = norm2(in);
|
||||||
|
RealD nr = norm2(resid);
|
||||||
|
|
||||||
|
std::cout<<GridLogMessage << "SchurRedBlackStag solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user