From a17684ebe2e30e0972f953a4996898b633ff2bc0 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Mon, 22 Jun 2015 12:49:44 +0100 Subject: [PATCH] Some small steps towards a multigrid --- lib/Make.inc | 2 +- lib/algorithms/CoarsenedMatrix.h | 6 +- lib/algorithms/iterative/ConjugateGradient.h | 2 +- .../PrecGeneralisedConjugateResidual.h | 6 +- lib/qcd/action/fermion/WilsonFermion5D.cc | 39 ++-- lib/qcd/utils/SUn.h | 4 +- tests/InvSqrt.gnu | 0 tests/Make.inc | 14 +- tests/Sqrt.gnu | 2 - tests/Test_cayley_ldop_cg.cc | 155 --------------- tests/Test_cayley_ldop_cr.cc | 103 ++-------- tests/Test_dwf_hdcr.cc | 188 ++++++++++++++++++ 12 files changed, 243 insertions(+), 278 deletions(-) delete mode 100644 tests/InvSqrt.gnu delete mode 100644 tests/Sqrt.gnu delete mode 100644 tests/Test_cayley_ldop_cg.cc create mode 100644 tests/Test_dwf_hdcr.cc diff --git a/lib/Make.inc b/lib/Make.inc index 6b5ae64f..91a06fcb 100644 --- a/lib/Make.inc +++ b/lib/Make.inc @@ -1,4 +1,4 @@ -HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h +HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/gauge/GaugeActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h index 2576a696..9ea77b84 100644 --- a/lib/algorithms/CoarsenedMatrix.h +++ b/lib/algorithms/CoarsenedMatrix.h @@ -85,10 +85,6 @@ namespace Grid { void Orthogonalise(void){ CoarseScalar InnerProd(CoarseGrid); blockOrthogonalise(InnerProd,subspace); -#if 1 - // CheckOrthogonal(); -#endif - } void CheckOrthogonal(void){ CoarseVector iProj(CoarseGrid); @@ -125,7 +121,7 @@ namespace Grid { RealD scale; - ConjugateGradient CG(1.0e-4,10000); + ConjugateGradient CG(1.0e-3,10000); FineField noise(FineGrid); FineField Mn(FineGrid); diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h index e281d765..3829dad8 100644 --- a/lib/algorithms/iterative/ConjugateGradient.h +++ b/lib/algorithms/iterative/ConjugateGradient.h @@ -15,7 +15,7 @@ public: Integer MaxIterations; int verbose; ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { - verbose=0; + verbose=1; }; diff --git a/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h b/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h index bd5bfbfe..77b7df49 100644 --- a/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h +++ b/lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h @@ -48,11 +48,11 @@ namespace Grid { if(cp=0)&&(dir<4) ); //must do x,y,z or t; + // assert( (disp==1)||(disp==-1) ); + // assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t; WilsonCompressor compressor(DaggerNo); Stencil.HaloExchange(in,comm_buf,compressor); @@ -100,7 +100,7 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int assert(dirdisp<=7); assert(dirdisp>=0); -PARALLEL_FOR_LOOP +//PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ for(int s=0;soSites();ss++){ - for(int s=0;soSites();ss++){ - for(int s=0;soSites();ss++){ for(int s=0;soSites();ss++){ for(int s=0;soSites();ss++){ + for(int ss=0;ssoSites();ss++){ subgroup._odata[ss]()()(0,0) = source._odata[ss]()()(i0,i0); subgroup._odata[ss]()()(0,1) = source._odata[ss]()()(i0,i1); subgroup._odata[ss]()()(1,0) = source._odata[ss]()()(i1,i0); @@ -201,7 +201,7 @@ PARALLEL_FOR_LOOP dest = 1.0; // start out with identity PARALLEL_FOR_LOOP - for(int ss=0;ss!=grid->oSites();ss++){ + for(int ss=0;ssoSites();ss++){ dest._odata[ss]()()(i0,i0) = subgroup._odata[ss]()()(0,0); dest._odata[ss]()()(i0,i1) = subgroup._odata[ss]()()(0,1); dest._odata[ss]()()(i1,i0) = subgroup._odata[ss]()()(1,0); diff --git a/tests/InvSqrt.gnu b/tests/InvSqrt.gnu deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/Make.inc b/tests/Make.inc index 4d636cc5..4ad511bd 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,15 +1,11 @@ -bin_PROGRAMS = Test_GaugeAction Test_Metropolis Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd +bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_GaugeAction_SOURCES=Test_GaugeAction.cc Test_GaugeAction_LDADD=-lGrid -Test_Metropolis_SOURCES=Test_Metropolis.cc -Test_Metropolis_LDADD=-lGrid - - Test_cayley_cg_SOURCES=Test_cayley_cg.cc Test_cayley_cg_LDADD=-lGrid @@ -22,10 +18,6 @@ Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc Test_cayley_even_odd_LDADD=-lGrid -Test_cayley_ldop_cg_SOURCES=Test_cayley_ldop_cg.cc -Test_cayley_ldop_cg_LDADD=-lGrid - - Test_cayley_ldop_cr_SOURCES=Test_cayley_ldop_cr.cc Test_cayley_ldop_cr_LDADD=-lGrid @@ -78,6 +70,10 @@ Test_dwf_fpgcr_SOURCES=Test_dwf_fpgcr.cc Test_dwf_fpgcr_LDADD=-lGrid +Test_dwf_hdcr_SOURCES=Test_dwf_hdcr.cc +Test_dwf_hdcr_LDADD=-lGrid + + Test_gamma_SOURCES=Test_gamma.cc Test_gamma_LDADD=-lGrid diff --git a/tests/Sqrt.gnu b/tests/Sqrt.gnu deleted file mode 100644 index ae56ab97..00000000 --- a/tests/Sqrt.gnu +++ /dev/null @@ -1,2 +0,0 @@ -f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866)); -f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422)); diff --git a/tests/Test_cayley_ldop_cg.cc b/tests/Test_cayley_ldop_cg.cc deleted file mode 100644 index f960ab7b..00000000 --- a/tests/Test_cayley_ldop_cg.cc +++ /dev/null @@ -1,155 +0,0 @@ -#include -#include -#include - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - - -template -struct scal { - d internal; -}; - - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; - -double lowpass(double x) -{ - return pow(x*x+1.0,-2); -} - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - Chebyshev filter(-150.0, 150.0,16, lowpass); - ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc); - filter.csv(csv); - csv.close(); - - const int Ls=8; - - GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi()); - GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); - - GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - - // Construct a coarsened grid - std::vector 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); - - LatticeFermion src(FGrid); gaussian(RNG5,src); - LatticeFermion result(FGrid); result=zero; - LatticeFermion ref(FGrid); ref=zero; - LatticeFermion tmp(FGrid); - LatticeFermion err(FGrid); - LatticeGaugeField Umu(UGrid); - - - - //gaussian(RNG4,Umu); - //random(RNG4,Umu); - - NerscField header; - std::string file("./ckpoint_lat.400"); - readNerscConfiguration(Umu,header,file); - // SU3::ColdConfiguration(RNG4,Umu); - // SU3::TepidConfiguration(RNG4,Umu); - // SU3::HotConfiguration(RNG4,Umu); - // Umu=zero; - -#if 0 - LatticeColourMatrix U(UGrid); - for(int nn=0;nn(Umu,nn); - U=U*adj(U)-1.0; - std::cout<<"SU3 test "< HermIndefOp(Ddwf); - - const int nbasis = 8; - -#if 0 - std::vector subspace(nbasis,FGrid); - LatticeFermion noise(FGrid); - LatticeFermion ms(FGrid); - for(int b=0;b HermDefOp(Ddwf); - ConjugateGradient CG(1.0e-4,10000); - - for(int i=0;i<1;i++){ - - CG(HermDefOp,noise,subspace[b]); - noise = subspace[b]; - - scale = pow(norm2(noise),-0.5); - noise=noise*scale; - HermDefOp.Op(noise,ms); std::cout << "filt "< "< HermDefOp(Ddwf); - typedef Aggregation Subspace; - Subspace Aggregates(Coarse5d,FGrid); - Aggregates.CreateSubspace(RNG5,HermDefOp); - std::cout << "Called aggregation class"<< std::endl; -#endif - - - typedef CoarsenedMatrix LittleDiracOperator; - typedef LittleDiracOperator::CoarseVector CoarseVector; - - LittleDiracOperator LittleDiracOp(*Coarse5d); - LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); - - CoarseVector c_src (Coarse5d); - CoarseVector c_res (Coarse5d); - gaussian(CRNG,c_src); - c_res=zero; - - std::cout << "Solving CG on coarse space "<< std::endl; - MdagMLinearOperator PosdefLdop(LittleDiracOp); - ConjugateGradient CG(1.0e-6,10000); - CG(PosdefLdop,c_src,c_res); - - std::cout << "Done "<< std::endl; - Grid_finalize(); -} diff --git a/tests/Test_cayley_ldop_cr.cc b/tests/Test_cayley_ldop_cr.cc index c4001b55..8926c19c 100644 --- a/tests/Test_cayley_ldop_cr.cc +++ b/tests/Test_cayley_ldop_cr.cc @@ -1,38 +1,13 @@ #include -#include -#include using namespace std; using namespace Grid; using namespace Grid::QCD; - -template -struct scal { - d internal; -}; - - Gamma::GammaMatrix Gmu [] = { - Gamma::GammaX, - Gamma::GammaY, - Gamma::GammaZ, - Gamma::GammaT - }; - -double lowpass(double x) -{ - return pow(x*x+1.0,-2); -} - int main (int argc, char ** argv) { Grid_init(&argc,&argv); - Chebyshev filter(-150.0, 150.0,16, lowpass); - ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc); - filter.csv(csv); - csv.close(); - const int Ls=8; GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi()); @@ -41,7 +16,9 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); - // Construct a coarsened grid + /////////////////////////////////////////////////// + // Construct a coarsened grid; utility for this? + /////////////////////////////////////////////////// std::vector clatt = GridDefaultLatt(); for(int d=0;d(Umu,nn); - U=U*adj(U)-1.0; - std::cout<<"SU3 test "< HermIndefOp(Ddwf); const int nbasis = 8; -#if 0 - std::vector subspace(nbasis,FGrid); - LatticeFermion noise(FGrid); - LatticeFermion ms(FGrid); - for(int b=0;b Subspace; + typedef CoarsenedMatrix LittleDiracOperator; + typedef LittleDiracOperator::CoarseVector CoarseVector; - gaussian(RNG5,noise); - RealD scale = pow(norm2(noise),-0.5); - noise=noise*scale; - - HermIndefOp.Op(noise,ms); std::cout << "Noise "< HermDefOp(Ddwf); - ConjugateGradient CG(1.0e-4,10000); - - for(int i=0;i<1;i++){ - - CG(HermDefOp,noise,subspace[b]); - noise = subspace[b]; - - scale = pow(norm2(noise),-0.5); - noise=noise*scale; - HermDefOp.Op(noise,ms); std::cout << "filt "< "< HermDefOp(Ddwf); - typedef Aggregation Subspace; Subspace Aggregates(Coarse5d,FGrid); Aggregates.CreateSubspace(RNG5,HermDefOp); - std::cout << "Called aggregation class"<< std::endl; -#endif - typedef CoarsenedMatrix LittleDiracOperator; - typedef LittleDiracOperator::CoarseVector CoarseVector; - LittleDiracOperator LittleDiracOp(*Coarse5d); LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); @@ -145,18 +80,22 @@ int main (int argc, char ** argv) gaussian(CRNG,c_src); c_res=zero; - - std::cout << "Solving CG on coarse space "<< std::endl; + std::cout << "**************************************************"<< std::endl; + std::cout << "Solving mdagm-CG on coarse space "<< std::endl; + std::cout << "**************************************************"<< std::endl; MdagMLinearOperator PosdefLdop(LittleDiracOp); ConjugateGradient CG(1.0e-6,10000); CG(PosdefLdop,c_src,c_res); - - std::cout << "Solving MCR on coarse space "<< std::endl; + std::cout << "**************************************************"<< std::endl; + std::cout << "Solving indef-MCR on coarse space "<< std::endl; + std::cout << "**************************************************"<< std::endl; HermitianLinearOperator HermIndefLdop(LittleDiracOp); ConjugateResidual MCR(1.0e-6,10000); MCR(HermIndefLdop,c_src,c_res); + std::cout << "**************************************************"<< std::endl; std::cout << "Done "<< std::endl; + std::cout << "**************************************************"<< std::endl; Grid_finalize(); } diff --git a/tests/Test_dwf_hdcr.cc b/tests/Test_dwf_hdcr.cc new file mode 100644 index 00000000..44b05714 --- /dev/null +++ b/tests/Test_dwf_hdcr.cc @@ -0,0 +1,188 @@ +#include +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +class MultiGridPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef CoarsenedMatrix CoarseOperator; + + typedef typename Aggregation::siteVector siteVector; + typedef typename Aggregation::CoarseScalar CoarseScalar; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + + Aggregates & _Aggregates; + CoarseOperator & _CoarseOperator; + FineOperator & _FineOperator; + + // Constructor + MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, FineOperator &Fine) + : _Aggregates(Agg), + _CoarseOperator(Coarse), + _FineOperator(Fine) + { + } + + void operator()(const FineField &in, FineField & out) { + + FineField Min(in._grid); + + CoarseVector Csrc(_CoarseOperator.Grid()); + CoarseVector Ctmp(_CoarseOperator.Grid()); + CoarseVector Csol(_CoarseOperator.Grid()); + + // Monitor completeness of low mode space + _Aggregates.ProjectToSubspace (Csrc,in); + _Aggregates.PromoteFromSubspace(Csrc,out); + std::cout<<"Completeness: "< MCR(1.0e-2,1000); + ConjugateGradient CG(1.0e-2,10000); + + // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] + // Smoothing step, followed by coarse grid correction + + MCR(_FineOperator,in,Min); + _FineOperator.Op(Min,out); + out = in -out; // out = in - A Min + + MdagMLinearOperator MdagMOp(_CoarseOperator); + HermitianLinearOperator HermOp(_CoarseOperator); + Csol=zero; + _Aggregates.ProjectToSubspace (Csrc,out); + HermOp.AdjOp(Csrc,Ctmp);// Normal equations + CG(MdagMOp ,Ctmp,Csol); + _Aggregates.PromoteFromSubspace(Csol,out); + + out = Min + out;; + } + +}; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + const int Ls=8; + + GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi()); + GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid); + + GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); + GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + + /////////////////////////////////////////////////// + // Construct a coarsened grid; utility for this? + /////////////////////////////////////////////////// + std::vector 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); + + LatticeFermion src(FGrid); gaussian(RNG5,src); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); + LatticeFermion err(FGrid); + LatticeGaugeField Umu(UGrid); + + NerscField header; + std::string file("./ckpoint_lat.4000"); + readNerscConfiguration(Umu,header,file); + + // SU3::ColdConfiguration(RNG4,Umu); + // SU3::TepidConfiguration(RNG4,Umu); + // SU3::HotConfiguration(RNG4,Umu); + // Umu=zero; + + RealD mass=0.04; + RealD M5=1.8; + + std::cout << "**************************************************"<< std::endl; + std::cout << "Building g5R5 hermitian DWF operator" < Subspace; + typedef CoarsenedMatrix CoarseOperator; + typedef CoarseOperator::CoarseVector CoarseVector; + + std::cout << "**************************************************"<< std::endl; + std::cout << "Calling Aggregation class to build subspace" < HermDefOp(Ddwf); + Subspace Aggregates(Coarse5d,FGrid); + Aggregates.CreateSubspace(RNG5,HermDefOp); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Building coarse representation of Indef operator" < HermIndefOp(Ddwf); + CoarsenedMatrix LDOp(*Coarse5d); + LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Testing some coarse space solvers " < PosdefLdop(LDOp); + ConjugateGradient CG(1.0e-6,10000); + CG(PosdefLdop,c_src,c_res); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Solving indef-MCR on coarse space "<< std::endl; + std::cout << "**************************************************"<< std::endl; + HermitianLinearOperator HermIndefLdop(LDOp); + ConjugateResidual MCR(1.0e-6,10000); + //MCR(HermIndefLdop,c_src,c_res); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Building deflation preconditioner "<< std::endl; + std::cout << "**************************************************"<< std::endl; + + MultiGridPreconditioner Precon(Aggregates, LDOp,HermIndefOp); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Building a one level PGCR "<< std::endl; + std::cout << "**************************************************"<< std::endl; + TrivialPrecon simple; + PrecGeneralisedConjugateResidual GCR(1.0e-6,10000,simple,8,64); + GCR(HermIndefOp,src,result); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Building a two level PGCR "<< std::endl; + std::cout << "**************************************************"<< std::endl; + PrecGeneralisedConjugateResidual PGCR(1.0e-6,10000,Precon,8,64); + PGCR(HermIndefOp,src,result); + + std::cout << "**************************************************"<< std::endl; + std::cout << "Done "<< std::endl; + std::cout << "**************************************************"<< std::endl; + Grid_finalize(); +}