portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-10 19:36:56 +01:00
Code Releases Activity

Made checkerboard choice in staggered preconditioned solvers switchable

Browse Source
This commit is contained in:
Chulwoo Jung
2018-02-13 11:49:02 -05:00
parent dd8cd8a8e8
commit 4ac85b3e8f
1 changed files with 38 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -175,11 +175,8 @@ namespace Grid {
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick // Wrap the usual normal equations Schur trick
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
SchurRedBlackStagSolve(OperatorFunction<Field> &HermitianRBSolver) : SchurRedBlackStagSolve(OperatorFunction<Field> &HermitianRBSolver, int cb) :
_HermitianRBSolver(HermitianRBSolver) _HermitianRBSolver(HermitianRBSolver), CBfactorise(cb) {}
{
CBfactorise=0;
};
template<class Matrix> template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){ void operator() (Matrix & _Matrix,const Field &in, Field &out){
@ -190,6 +187,8 @@ namespace Grid {
GridBase *fgrid= _Matrix.Grid(); GridBase *fgrid= _Matrix.Grid();
SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix); SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix);
int Schur = CBfactorise;
int Other = 1 - CBfactorise;
Field src_e(grid); Field src_e(grid);
Field src_o(grid); Field src_o(grid);
@ -199,39 +198,39 @@ namespace Grid {
Field Mtmp(grid); Field Mtmp(grid);
Field resid(fgrid); Field resid(fgrid);
pickCheckerboard(Even,src_e,in); pickCheckerboard(Other,src_e,in);
pickCheckerboard(Odd ,src_o,in); pickCheckerboard(Schur ,src_o,in);
pickCheckerboard(Even,sol_e,out); pickCheckerboard(Other,sol_e,out);
pickCheckerboard(Odd ,sol_o,out); pickCheckerboard(Schur ,sol_o,out);
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = Mdag * (source_o - Moe MeeInv source_e) // src_o = Mdag * (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Other);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Schur);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.checkerboard ==Schur);
#if 0 #if 0
// get the right MpcDag // get the right MpcDag
// _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); // _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Schur);
#else #else
_Matrix.Mooee(tmp,src_o); assert(src_o.checkerboard ==Odd); _Matrix.Mooee(tmp,src_o); assert(src_o.checkerboard ==Schur);
#endif #endif
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Schur);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Other);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.checkerboard ==Other);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Other);
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Other);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Schur );
// Verify the unprec residual // Verify the unprec residual
_Matrix.M(out,resid); _Matrix.M(out,resid);
@ -489,11 +488,8 @@ namespace Grid {
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// Wrap the usual normal equations Schur trick // Wrap the usual normal equations Schur trick
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
SchurRedBlackStagMixed(LinearFunction<Field> &HermitianRBSolver) : SchurRedBlackStagMixed(LinearFunction<Field> &HermitianRBSolver, int cb) :
_HermitianRBSolver(HermitianRBSolver) _HermitianRBSolver(HermitianRBSolver), CBfactorise(cb) {}
{
CBfactorise=0;
};
template<class Matrix> template<class Matrix>
void operator() (Matrix & _Matrix,const Field &in, Field &out){ void operator() (Matrix & _Matrix,const Field &in, Field &out){
@ -504,6 +500,8 @@ namespace Grid {
GridBase *fgrid= _Matrix.Grid(); GridBase *fgrid= _Matrix.Grid();
SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix); SchurStagOperator<Matrix,Field> _HermOpEO(_Matrix);
int Schur = CBfactorise;
int Other = 1 - CBfactorise;
Field src_e(grid); Field src_e(grid);
Field src_o(grid); Field src_o(grid);
@ -513,35 +511,35 @@ namespace Grid {
Field Mtmp(grid); Field Mtmp(grid);
Field resid(fgrid); Field resid(fgrid);
pickCheckerboard(Even,src_e,in); pickCheckerboard(Other,src_e,in);
pickCheckerboard(Odd ,src_o,in); pickCheckerboard(Schur ,src_o,in);
pickCheckerboard(Even,sol_e,out); pickCheckerboard(Other,sol_e,out);
pickCheckerboard(Odd ,sol_o,out); pickCheckerboard(Schur ,sol_o,out);
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
// src_o = Mdag * (source_o - Moe MeeInv source_e) // src_o = Mdag * (source_o - Moe MeeInv source_e)
///////////////////////////////////////////////////// /////////////////////////////////////////////////////
_Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Even); _Matrix.MooeeInv(src_e,tmp); assert( tmp.checkerboard ==Other);
_Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Odd); _Matrix.Meooe (tmp,Mtmp); assert( Mtmp.checkerboard ==Schur);
tmp=src_o-Mtmp; assert( tmp.checkerboard ==Odd); tmp=src_o-Mtmp; assert( tmp.checkerboard ==Schur);
_Matrix.Mooee(tmp,src_o); assert(src_o.checkerboard ==Odd); _Matrix.Mooee(tmp,src_o); assert(src_o.checkerboard ==Schur);
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
// Call the red-black solver // Call the red-black solver
////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl; std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
// _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd); // _HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Schur);
_HermitianRBSolver(src_o,sol_o); assert(sol_o.checkerboard==Odd); _HermitianRBSolver(src_o,sol_o); assert(sol_o.checkerboard==Other);
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
// sol_e = M_ee^-1 * ( src_e - Meo sol_o )... // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
/////////////////////////////////////////////////// ///////////////////////////////////////////////////
_Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Even); _Matrix.Meooe(sol_o,tmp); assert( tmp.checkerboard ==Other);
src_e = src_e-tmp; assert( src_e.checkerboard ==Even); src_e = src_e-tmp; assert( src_e.checkerboard ==Other);
_Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Even); _Matrix.MooeeInv(src_e,sol_e); assert( sol_e.checkerboard ==Other);
setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Even); setCheckerboard(out,sol_e); assert( sol_e.checkerboard ==Other);
setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Odd ); setCheckerboard(out,sol_o); assert( sol_o.checkerboard ==Schur );
// Verify the unprec residual // Verify the unprec residual
_Matrix.M(out,resid); _Matrix.M(out,resid);