From b9dcad89e8508c026421acde95de9e1eb190a956 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 7 Sep 2023 10:53:22 -0400 Subject: [PATCH] Test cases for coarsening with non-local stencil --- tests/debug/Test_general_coarse.cc | 15 +- tests/debug/Test_general_coarse_hdcg.cc | 203 ++++++++++++++++++++++ tests/debug/Test_general_coarse_pvdagm.cc | 66 +++++-- 3 files changed, 263 insertions(+), 21 deletions(-) create mode 100644 tests/debug/Test_general_coarse_hdcg.cc diff --git a/tests/debug/Test_general_coarse.cc b/tests/debug/Test_general_coarse.cc index 169f9378..1430ed1d 100644 --- a/tests/debug/Test_general_coarse.cc +++ b/tests/debug/Test_general_coarse.cc @@ -81,14 +81,11 @@ public: } void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); } void HermOp(const Field &in, Field &out){ - std::cout << "HermOp"< + + 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 +#include +#include +#include + +using namespace std; +using namespace Grid; + +// Want Op in CoarsenOp to call MatPcDagMatPc +template +class HermOpAdaptor : public LinearOperatorBase +{ + LinearOperatorBase & wrapped; +public: + HermOpAdaptor(LinearOperatorBase &wrapme) : wrapped(wrapme) {}; + void OpDiag (const Field &in, Field &out) { assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); }; + void Op (const Field &in, Field &out){ + wrapped.HermOp(in,out); + } + void AdjOp (const Field &in, Field &out){ + wrapped.HermOp(in,out); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); } + void HermOp(const Field &in, Field &out){ + wrapped.HermOp(in,out); + } + +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=16; + + 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); + + // Construct a coarsened grid + // 4^4 cell + Coordinate clatt = GridDefaultLatt(); + for(int d=0;d seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + std::vector cseeds({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); + + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + std::string file("ckpoint_lat.4000"); + NerscIO::readConfiguration(Umu,header,file); + + RealD mass=0.01; + RealD M5=1.8; + + RealD b=1.5; + RealD c=0.5; + MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + MobiusFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5,b,c); + + const int nbasis = 4; + const int cb = 0 ; + LatticeFermion prom(FrbGrid); + + typedef GeneralCoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; + + NextToNextToNextToNearestStencilGeometry5D geom; + + std::cout< HermOpEO(Ddwf); + HermOpAdaptor HOA(HermOpEO); + + // Run power method on HOA?? + LatticeFermion result(FrbGrid); result=Zero(); + LatticeFermion ref(FrbGrid); ref=Zero(); + LatticeFermion tmp(FrbGrid); + LatticeFermion err(FrbGrid); + + { + LatticeFermion src(FrbGrid); random(RNG5,src); + PowerMethod PM; PM(HermOpEO,src); + } + // exit(0); + + // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp + typedef Aggregation Subspace; + Subspace Aggregates(Coarse5d,FrbGrid,cb); + Aggregates.CreateSubspaceChebyshev(RNG5, + HermOpEO, + nbasis, + 90.0, + 0.1, + 500, + 500, + 100, + 0.0); + //////////////////////////////////////////////////////////// + // Need to check about red-black grid coarsening + //////////////////////////////////////////////////////////// + LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d); + LittleDiracOp.CoarsenOperator(HOA,Aggregates); + + std::cout< subspace(nbasis,FrbGrid); + subspace=Aggregates.subspace; + + Complex one(1.0); + c_src = one; // 1 in every element for vector 1. + blockPromote(c_src,err,subspace); + + prom=Zero(); + for(int b=0;b Guess; + RealD tol = 1.0e-8; + int maxit=2000; + ConjugateGradient CG(tol,maxit,false); + HermitianLinearOperator Hop (LittleDiracOp); + CG(Hop, c_src, c_res); + Grid_finalize(); + return 0; +} diff --git a/tests/debug/Test_general_coarse_pvdagm.cc b/tests/debug/Test_general_coarse_pvdagm.cc index c9652a2c..f834b3bd 100644 --- a/tests/debug/Test_general_coarse_pvdagm.cc +++ b/tests/debug/Test_general_coarse_pvdagm.cc @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -141,7 +141,7 @@ int main (int argc, char ** argv) { Grid_init(&argc,&argv); - const int Ls=16; + const int Ls=2; GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); @@ -152,7 +152,7 @@ int main (int argc, char ** argv) // Construct a coarsened grid Coordinate clatt = GridDefaultLatt(); for(int d=0;d Subspace; Subspace AggregatesPD(Coarse5d,FGrid,cb); AggregatesPD.CreateSubspaceChebyshev(RNG5, - HOA, - nbasis, - 5000.0, - 0.02, - 100, - 50, - 50, - 0.0); + HOA, + nbasis, + 5000.0, + 0.02, + 100, + 50, + 50, + 0.0); LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d); LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD); + std::cout< subspace(nbasis,FGrid); + subspace=AggregatesPD.subspace; + + Complex one(1.0); + c_src = one; // 1 in every element for vector 1. + blockPromote(c_src,err,subspace); + + prom=Zero(); + for(int b=0;b