From 1e54882f7145bd38db8ac1681cb7d4f9bceb2297 Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Wed, 4 Oct 2017 10:51:06 +0100 Subject: [PATCH 1/3] Stagger --- tests/solver/Test_staggered_cg_prec.cc | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/solver/Test_staggered_cg_prec.cc b/tests/solver/Test_staggered_cg_prec.cc index 66f11d3d..0a803c21 100644 --- a/tests/solver/Test_staggered_cg_prec.cc +++ b/tests/solver/Test_staggered_cg_prec.cc @@ -83,5 +83,14 @@ int main (int argc, char ** argv) ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,res_o); + FermionField tmp(&RBGrid); + + HermOpEO.Mpc(res_o,tmp); + std::cout << "check Mpc resid " << axpy_norm(tmp,-1.0,src_o,tmp)/norm2(src_o) << "\n"; + + RealD n1,n2; + HermOpEO.MpcDagMpc(res_o,tmp,n1,n2); + std::cout << "check MpcDagMpc resid " << axpy_norm(tmp,-1.0,src_o,tmp)/norm2(src_o) << "\n"; + Grid_finalize(); } From f0e084a88c98c2ebf94a7c5f9f6f145699bfd827 Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Tue, 10 Oct 2017 10:00:43 +0100 Subject: [PATCH 2/3] Schur staggered --- lib/algorithms/LinearOperator.h | 104 +++++++++- lib/algorithms/iterative/SchurRedBlack.h | 240 +++++++++++++++++++++++ 2 files changed, 342 insertions(+), 2 deletions(-) diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index 6cb77296..6e4da248 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -192,10 +192,10 @@ namespace Grid { ni=Mpc(in,tmp); no=MpcDag(tmp,out); } - void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ MpcDagMpc(in,out,n1,n2); } - void HermOp(const Field &in, Field &out){ + virtual void HermOp(const Field &in, Field &out){ RealD n1,n2; HermOpAndNorm(in,out,n1,n2); } @@ -300,6 +300,106 @@ namespace Grid { } }; + // + template + class SchurStaggeredOperator : public SchurOperatorBase { + protected: + Matrix &_Mat; + public: + SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){}; + virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + ComplexD dot; + n2=Mpc(in,out); + dot= innerProduct(in,out); + n1= real(dot); + } + virtual void HermOp(const Field &in, Field &out){ + Mpc(in,out); + } + virtual RealD Mpc (const Field &in, Field &out) { + Field tmp(in._grid); + _Mat.Meooe(in,tmp); + _Mat.MooeeInv(tmp,out); + _Mat.MeooeDag(out,tmp); + _Mat.Mooee(in,out); + return axpy_norm(out,-1.0,tmp,out); + } + virtual RealD MpcDag (const Field &in, Field &out){ + return Mpc(in,out); + } + virtual void MpcDagMpc(const Field &in, Field &out,RealD &ni,RealD &no) { + assert(0);// Never need with staggered + } + }; + template using SchurStagOperator = SchurStaggeredOperator; + + // This is specific to (Z)mobius fermions + template + class KappaSimilarityTransform { + public: + + typedef typename Matrix::Coeff_t Coeff_t; + std::vector kappa, kappaDag, kappaInv, kappaInvDag; + + KappaSimilarityTransform (Matrix &zmob) { + for (int i=0;i<(int)zmob.bs.size();i++) { + Coeff_t k = 1.0 / ( 2.0 * (zmob.bs[i] *(4 - zmob.M5) + 1.0) ); + kappa.push_back( k ); + kappaDag.push_back( conj(k) ); + kappaInv.push_back( 1.0 / k ); + kappaInvDag.push_back( 1.0 / conj(k) ); + } + } + + template + void sscale(const Lattice& in, Lattice& out, Coeff_t* s) { + GridBase *grid=out._grid; + out.checkerboard = in.checkerboard; + assert(grid->_simd_layout[0] == 1); // should be fine for ZMobius for now + int Ls = grid->_rdimensions[0]; + parallel_for(int ss=0;ssoSites();ss++){ + vobj tmp = s[ss % Ls]*in._odata[ss]; + vstream(out._odata[ss],tmp); + } + } + + RealD sscale_norm(const Field& in, Field& out, Coeff_t* s) { + sscale(in,out,s); + return norm2(out); + } + + virtual RealD M (const Field& in, Field& out) { return sscale_norm(in,out,&kappa[0]); } + virtual RealD MDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaDag[0]);} + virtual RealD MInv (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInv[0]);} + virtual RealD MInvDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInvDag[0]);} + + }; + + template + class SchurDiagTwoKappaOperator : public SchurOperatorBase { + public: + KappaSimilarityTransform _S; + SchurDiagTwoOperator _Mat; + + SchurDiagTwoKappaOperator (Matrix &Mat): _S(Mat), _Mat(Mat) {}; + + virtual RealD Mpc (const Field &in, Field &out) { + Field tmp(in._grid); + + _S.MInv(in,out); + _Mat.Mpc(out,tmp); + return _S.M(tmp,out); + + } + virtual RealD MpcDag (const Field &in, Field &out){ + Field tmp(in._grid); + + _S.MDag(in,out); + _Mat.MpcDag(out,tmp); + return _S.MInvDag(tmp,out); + } + }; + ///////////////////////////////////////////////////////////// // Base classes for functions of operators diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h index 5caabb4b..b6eab762 100644 --- a/lib/algorithms/iterative/SchurRedBlack.h +++ b/lib/algorithms/iterative/SchurRedBlack.h @@ -63,6 +63,85 @@ Author: Peter Boyle */ 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 SchurRedBlackStaggeredSolve { + private: + OperatorFunction & _HermitianRBSolver; + int CBfactorise; + public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + SchurRedBlackStaggeredSolve(OperatorFunction &HermitianRBSolver) : + _HermitianRBSolver(HermitianRBSolver) + { + CBfactorise=0; + }; + + template + 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(); + + SchurStaggeredOperator _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< using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve; + /////////////////////////////////////////////////////////////////////////////////////////////////////// // Take a matrix and form a Red Black solver calling a Herm solver // Use of RB info prevents making SchurRedBlackSolve conform to standard interface @@ -141,5 +220,166 @@ 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 SchurRedBlackDiagTwoSolve { + private: + OperatorFunction & _HermitianRBSolver; + int CBfactorise; + public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + SchurRedBlackDiagTwoSolve(OperatorFunction &HermitianRBSolver) : + _HermitianRBSolver(HermitianRBSolver) + { + CBfactorise=0; + }; + + template + 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(); + + SchurDiagTwoOperator _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); + + // get the right MpcDag + _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + + ////////////////////////////////////////////////////////////// + // Call the red-black solver + ////////////////////////////////////////////////////////////// + std::cout< class SchurRedBlackDiagTwoMixed { + private: + LinearFunction & _HermitianRBSolver; + int CBfactorise; + public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + SchurRedBlackDiagTwoMixed(LinearFunction &HermitianRBSolver) : + _HermitianRBSolver(HermitianRBSolver) + { + CBfactorise=0; + }; + + template + 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(); + + SchurDiagTwoOperator _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); + + // get the right MpcDag + _HermOpEO.MpcDag(tmp,src_o); assert(src_o.checkerboard ==Odd); + + ////////////////////////////////////////////////////////////// + // Call the red-black solver + ////////////////////////////////////////////////////////////// + std::cout< Date: Tue, 10 Oct 2017 12:02:18 +0100 Subject: [PATCH 3/3] Schur for staggered --- lib/algorithms/LinearOperator.h | 80 +----------- tests/solver/Test_staggered_block_cg_prec.cc | 130 +++++++++++++++++++ tests/solver/Test_staggered_cg_prec.cc | 6 +- 3 files changed, 135 insertions(+), 81 deletions(-) create mode 100644 tests/solver/Test_staggered_block_cg_prec.cc diff --git a/lib/algorithms/LinearOperator.h b/lib/algorithms/LinearOperator.h index 6e4da248..d402c5b7 100644 --- a/lib/algorithms/LinearOperator.h +++ b/lib/algorithms/LinearOperator.h @@ -162,15 +162,10 @@ namespace Grid { _Mat.M(in,out); } void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ - ComplexD dot; - _Mat.M(in,out); - dot= innerProduct(in,out); - n1=real(dot); - - dot = innerProduct(out,out); - n2=real(dot); + ComplexD dot= innerProduct(in,out); n1=real(dot); + n2=norm2(out); } void HermOp(const Field &in, Field &out){ _Mat.M(in,out); @@ -309,9 +304,9 @@ namespace Grid { SchurStaggeredOperator (Matrix &Mat): _Mat(Mat){}; virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ ComplexD dot; - n2=Mpc(in,out); + n2 = Mpc(in,out); dot= innerProduct(in,out); - n1= real(dot); + n1 = real(dot); } virtual void HermOp(const Field &in, Field &out){ Mpc(in,out); @@ -333,73 +328,6 @@ namespace Grid { }; template using SchurStagOperator = SchurStaggeredOperator; - // This is specific to (Z)mobius fermions - template - class KappaSimilarityTransform { - public: - - typedef typename Matrix::Coeff_t Coeff_t; - std::vector kappa, kappaDag, kappaInv, kappaInvDag; - - KappaSimilarityTransform (Matrix &zmob) { - for (int i=0;i<(int)zmob.bs.size();i++) { - Coeff_t k = 1.0 / ( 2.0 * (zmob.bs[i] *(4 - zmob.M5) + 1.0) ); - kappa.push_back( k ); - kappaDag.push_back( conj(k) ); - kappaInv.push_back( 1.0 / k ); - kappaInvDag.push_back( 1.0 / conj(k) ); - } - } - - template - void sscale(const Lattice& in, Lattice& out, Coeff_t* s) { - GridBase *grid=out._grid; - out.checkerboard = in.checkerboard; - assert(grid->_simd_layout[0] == 1); // should be fine for ZMobius for now - int Ls = grid->_rdimensions[0]; - parallel_for(int ss=0;ssoSites();ss++){ - vobj tmp = s[ss % Ls]*in._odata[ss]; - vstream(out._odata[ss],tmp); - } - } - - RealD sscale_norm(const Field& in, Field& out, Coeff_t* s) { - sscale(in,out,s); - return norm2(out); - } - - virtual RealD M (const Field& in, Field& out) { return sscale_norm(in,out,&kappa[0]); } - virtual RealD MDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaDag[0]);} - virtual RealD MInv (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInv[0]);} - virtual RealD MInvDag (const Field& in, Field& out) { return sscale_norm(in,out,&kappaInvDag[0]);} - - }; - - template - class SchurDiagTwoKappaOperator : public SchurOperatorBase { - public: - KappaSimilarityTransform _S; - SchurDiagTwoOperator _Mat; - - SchurDiagTwoKappaOperator (Matrix &Mat): _S(Mat), _Mat(Mat) {}; - - virtual RealD Mpc (const Field &in, Field &out) { - Field tmp(in._grid); - - _S.MInv(in,out); - _Mat.Mpc(out,tmp); - return _S.M(tmp,out); - - } - virtual RealD MpcDag (const Field &in, Field &out){ - Field tmp(in._grid); - - _S.MDag(in,out); - _Mat.MpcDag(out,tmp); - return _S.MInvDag(tmp,out); - } - }; - ///////////////////////////////////////////////////////////// // Base classes for functions of operators diff --git a/tests/solver/Test_staggered_block_cg_prec.cc b/tests/solver/Test_staggered_block_cg_prec.cc new file mode 100644 index 00000000..1d0117e0 --- /dev/null +++ b/tests/solver/Test_staggered_block_cg_prec.cc @@ -0,0 +1,130 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_unprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField; + typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField; + typename ImprovedStaggeredFermion5DR::ImplParams params; + + const int Ls=8; + + Grid_init(&argc,&argv); + + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(UGrid ); pRNG.SeedFixedIntegers(seeds); + GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds); + + FermionField src(FGrid); random(pRNG5,src); + FermionField src_o(FrbGrid); pickCheckerboard(Odd,src_o,src); + FermionField result_o(FrbGrid); result_o=zero; + RealD nrm = norm2(src); + + LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu); + + RealD mass=0.003; + ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass); + SchurDiagMooeeOperator HermOp(Ds); + + ConjugateGradient CG(1.0e-8,10000); + int blockDim = 0; + BlockConjugateGradient BCGrQ(BlockCGrQ,blockDim,1.0e-8,10000); + BlockConjugateGradient BCG (BlockCG,blockDim,1.0e-8,10000); + BlockConjugateGradient mCG (CGmultiRHS,blockDim,1.0e-8,10000); + + std::cout << GridLogMessage << "****************************************************************** "< HermOp4d(Ds4d); + FermionField src4d(UGrid); random(pRNG,src4d); + FermionField src4d_o(UrbGrid); pickCheckerboard(Odd,src4d_o,src4d); + FermionField result4d_o(UrbGrid); + + result4d_o=zero; + CG(HermOp4d,src4d_o,result4d_o); + std::cout << GridLogMessage << "************************************************************************ "<