/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_dwf_hdcr.cc Copyright (C) 2015 Author: Antonin Portelli Author: Peter Boyle Author: paboyle 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 using namespace std; using namespace Grid; class myclass: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(myclass, int, domaindecompose, int, domainsize, int, order, int, Ls, double, mq, double, lo, double, hi, int, steps); myclass(){}; }; RealD InverseApproximation(RealD x){ return 1.0/x; } 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; Matrix & _FineMatrix; FineOperator & _FineOperator; Matrix & _SmootherMatrix; FineOperator & _SmootherOperator; Guesser & _Guess; double cheby_hi; double cheby_lo; int cheby_ord; myclass _params; // Constructor MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, FineOperator &Fine,Matrix &FineMatrix, FineOperator &Smooth,Matrix &SmootherMatrix, Guesser &Guess_, myclass params_) : _Aggregates(Agg), _CoarseOperator(Coarse), _FineOperator(Fine), _FineMatrix(FineMatrix), _SmootherOperator(Smooth), _SmootherMatrix(SmootherMatrix), _Guess(Guess_), _params(params_) { } void PowerMethod(const FineField &in) { FineField p1(in.Grid()); FineField p2(in.Grid()); MdagMLinearOperator fMdagMOp(_FineMatrix); p1=in; for(int i=0;i<50;i++){ RealD absp1=std::sqrt(norm2(p1)); fMdagMOp.HermOp(p1,p2);// this is the G5 herm bit // _FineOperator.Op(p1,p2);// this is the G5 herm bit RealD absp2=std::sqrt(norm2(p2)); if(i%10==9) std::cout< CG(1.0e-3,100,false); ConjugateGradient fCG(1.0e-3,10,false); HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); MdagMLinearOperator fMdagMOp(_FineMatrix); FineField tmp(in.Grid()); FineField res(in.Grid()); FineField Min(in.Grid()); // Monitor completeness of low mode space _Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.PromoteFromSubspace(Csrc,out); std::cout< CG(1.0e-10,100000); ConjugateGradient fCG(1.0e-3,1000); HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); ShiftedMdagMLinearOperator fMdagMOp(_FineMatrix,0.1); FineField tmp(in.Grid()); FineField res(in.Grid()); FineField Qin(in.Grid()); // Monitor completeness of low mode space // _Aggregates.ProjectToSubspace (Csrc,in); // _Aggregates.PromoteFromSubspace(Csrc,out); // std::cout< fMdagMOp(_FineMatrix); ShiftedMdagMLinearOperator fMdagMOp(_SmootherMatrix,0.0); RealD Ni,r; Ni = norm2(in); for(int ilo=0;ilo<3;ilo++){ for(int ord=5;ord<50;ord*=2){ std::cout << " lo "< Cheby (lo[ilo],70.0,ord,InverseApproximation); Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M _FineOperator.Op(vec2,vec1);// this is the G5 herm bit vec1 = in - vec1; // tmp = in - A Min r=norm2(vec1); std::cout< CG(5.0e-2,100000); HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); // MdagMLinearOperator fMdagMOp(_FineMatrix); ShiftedMdagMLinearOperator fMdagMOp(_SmootherMatrix,0.0); FineField vec1(in.Grid()); FineField vec2(in.Grid()); Chebyshev Cheby (_params.lo,_params.hi,_params.order,InverseApproximation); Chebyshev ChebyAccu(_params.lo,_params.hi,_params.order,InverseApproximation); // _Aggregates.ProjectToSubspace (Csrc,in); // _Aggregates.PromoteFromSubspace(Csrc,out); // std::cout< block ({2,2,2,2}); const int nbasis= 32; auto 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); Gamma g5(Gamma::Algebra::Gamma5); LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; LatticeFermion result(FGrid); result=Zero(); LatticeFermion ref(FGrid); ref=Zero(); LatticeFermion tmp(FGrid); LatticeFermion err(FGrid); LatticeGaugeField Umu(UGrid); LatticeGaugeField UmuDD(UGrid); LatticeColourMatrix U(UGrid); LatticeColourMatrix zz(UGrid); FieldMetaData header; std::string file("./ckpoint_lat.4000"); NerscIO::readConfiguration(Umu,header,file); if ( params.domaindecompose ) { Lattice > coor(UGrid); zz=Zero(); for(int mu=0;mu(Umu,mu); U = where(mod(coor,params.domainsize)==(Integer)0,zz,U); PokeIndex(UmuDD,U,mu); } } else { UmuDD = Umu; } // SU3::ColdConfiguration(RNG4,Umu); // SU3::TepidConfiguration(RNG4,Umu); // SU3::HotConfiguration(RNG4,Umu); // Umu=Zero(); RealD mass=params.mq; RealD M5=1.8; std::cout< Subspace; typedef CoarsenedMatrix CoarseOperator; typedef CoarseOperator::CoarseVector CoarseVector; typedef CoarseOperator::siteVector siteVector; std::cout< HermDefOp(Ddwf); Subspace Aggregates(Coarse5d,FGrid,0); assert ( (nbasis & 0x1)==0); { int nb=nbasis/2; std::cout< HermIndefOp(Ddwf); Gamma5R5HermitianLinearOperator HermIndefOpDD(DdwfDD); CoarsenedMatrix LDOp(*Coarse5d,1); // Hermitian matrix LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); exit(0); CoarseVector c_src (Coarse5d); CoarseVector c_res (Coarse5d); gaussian(CRNG,c_src); result=Zero(); c_res=Zero(); ////////////////////////////////////////////////// // Deflate the course space. Recursive multigrid? ////////////////////////////////////////////////// typedef CoarsenedMatrix Level1Op; typedef CoarsenedMatrix,nbasis> Level2Op; auto cclatt = clatt; for(int d=0;d,nbasis> CoarseSubspace; CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); double c_first = 0.2; double c_div = 1.2; std::vector c_lo(nbasis/2); c_lo[0] = c_first; for(int b=1;b c_ord(nbasis/2,200); c_ord[0]=500; #define RECURSIVE_MULTIGRID #ifdef RECURSIVE_MULTIGRID std::cout< PosdefLdop(LDOp); // CoarseAggregates.CreateSubspaceChebyshev(CRNG,PosdefLdop,nbasis,14.0,c_lo,c_ord); // CoarseAggregates.CreateSubspaceRandom(CRNG); // Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix // HermitianLinearOperator L1LinOp(LDOp); // L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates); #endif // std::cout< simple; // ConjugateGradient fCG(1.0e-8,100000); // fCG(HermDefOp,src,result); std::cout< HermOpEO(Ddwf); ConjugateGradient pCG(1.0e-8,10000); // pCG(HermOpEO,src_o,result_o); std::cout< IRLHermOp(LDOp); Chebyshev IRLCheby(0.005,16.0,51); // IRLCheby.InitLowPass(0.01,18.0,51); FunctionHermOp IRLOpCheby(IRLCheby,IRLHermOp); PlainHermOp IRLOp (IRLHermOp); int Nstop=24; int Nk=24; int Nm=48; ImplicitlyRestartedLanczos IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-3,20); int Nconv; std::vector eval(Nm); std::vector evec(Nm,Coarse5d); IRL.calc(eval,evec,c_src,Nconv); std::cout< CG(3.0e-3,100000); // CG(PosdefLdop,c_src,c_res); std::cout< DeflCoarseGuesser(evec,eval); DeflCoarseGuesser(c_src,c_res); // CG(PosdefLdop,c_src,c_res); std::cout< CoarseZeroGuesser; MultiGridPreconditioner > Precon (Aggregates, LDOp, HermIndefOp,Ddwf, HermIndefOp,Ddwf, CoarseZeroGuesser, params); // Precon.PowerMethod(src); /* std::cout<,nbasis,Level1Op,ZeroGuesser > CoarsePrecon (CoarseAggregates, L2Op, L1LinOp,LDOp, L1LinOp,LDOp, CoarseZeroGuesser, cparams); CoarsePrecon.PowerMethod(c_src); */ /* std::cout< PGCR(1.0e-8,100000,Precon,8,8); std::cout< > DeflatedPrecon (Aggregates, LDOp, HermIndefOp,Ddwf, HermIndefOp,Ddwf, DeflCoarseGuesser, params); PrecGeneralisedConjugateResidual deflPGCR(1.0e-8,100000,DeflatedPrecon,16,16); std::cout< CPGCR(1.0e-3,10000,CoarsePrecon,8,8); std::cout<