diff --git a/tests/solver/Test_dwf_hdcr_16_rb.cc b/tests/solver/Test_dwf_hdcr_16_rb.cc new file mode 100644 index 00000000..b7900b04 --- /dev/null +++ b/tests/solver/Test_dwf_hdcr_16_rb.cc @@ -0,0 +1,397 @@ +/************************************************************************************* + + 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 + +using namespace std; +using namespace Grid; +/* Params + * Grid: + * block1(4) + * block2(4) + * + * Subspace + * * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100 + * * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100 + + * Smoother: + * * Fine: Cheby(hi, lo, order) -- 60,0.5,10 + * * Coarse: Cheby(hi, lo, order) -- 12,0.1,4 + + * Lanczos: + * CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61 + */ +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class SolverWrapper : public LinearFunction { +private: + CheckerBoardedSparseMatrixBase & _Matrix; + SchurRedBlackBase & _Solver; +public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations trick + ///////////////////////////////////////////////////// + SolverWrapper(CheckerBoardedSparseMatrixBase &Matrix, + SchurRedBlackBase &Solver) + : _Matrix(Matrix), _Solver(Solver) {}; + + void operator() (const Field &in, Field &out){ + + _Solver(_Matrix,in,out); // Mdag M out = Mdag in + + } +}; + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template class MirsSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & SmootherMatrix; + FineOperator & SmootherOperator; + RealD tol; + RealD shift; + int maxit; + + MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : + shift(_shift),tol(_tol),maxit(_maxit), + SmootherOperator(_SmootherOperator), + SmootherMatrix(_SmootherMatrix) + {}; + + void operator() (const Field &in, Field &out) + { + ZeroGuesser Guess; + ConjugateGradient CG(tol,maxit,false); + + Field src(in.Grid()); + + ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); + SmootherOperator.AdjOp(in,src); + Guess(src,out); + CG(MdagMOp,src,out); + } +}; + +template +class MultiGridPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef CoarsenedMatrix CoarseOperator; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + + Aggregates & _Aggregates; + CoarseOperator & _CoarseOperator; + Matrix & _FineMatrix; + FineOperator & _FineOperator; + Guesser & _Guess; + FineSmoother & _Smoother; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + +#define GridLogLevel std::cout << GridLogMessage < block ({2,2,2,2}); + std::vector blockc ({2,2,2,2}); + const int nbasis= 32; + const int nbasisc= 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); + LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; + LatticeFermion result(FGrid); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + std::string file("./ckpoint_lat.4000"); + //std::string file("./ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + 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; + LatticeFermion A(FGrid); + LatticeFermion B(FGrid); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0); + Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,200,150,0.0);// + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster + + for(int n=0;n Level1Op; + typedef CoarsenedMatrix,nbasisc> Level2Op; + + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); + + + GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d); + std::cout << " Making 5D coarse RB grid " <,nbasisc> CoarseSubspace; + // CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< PosdefLdop(LDOp); + typedef Level2Op::CoarseVector CoarseCoarseVector; + CoarseVector c_src(Coarse5d); c_src=1.0; + + std::cout< , SolverWrapper > TwoLevelMG; + typedef MultiGridPreconditioner,nbasisc,Level1Op, DeflatedGuesser, NormalEquations > CoarseMG; + typedef MultiGridPreconditioner, LinearFunction > ThreeLevelMG; + + ChebyshevSmoother FineSmoother(0.5,60.0,12,HermIndefOp,Ddwf); + std::cout< CoarseZeroGuesser; + ConjugateGradient CoarseCG(0.005,1000); + // SchurDiagMooeeOperator CoarseMpcDagMpc(LDOp); + SchurRedBlackDiagMooeeSolve CoarseRBCG(CoarseCG); + SolverWrapper CoarseSolver(LDOp,CoarseRBCG); + + // NormalEquations CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser); + TwoLevelMG TwoLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + CoarseSolver); + TwoLevelPrecon.Level(1); + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16); + l1PGCR.Level(1); + l1PGCR(src,result); + + std::cout< pCG(1.0e-8,60000); + result=Zero(); + // pCG(HermDefOp,src,result); + + std::cout< HermOpEO(Ddwf); + // pCG(HermOpEO,src_o,result_o); + + std::cout< PM; PM(HermDefOp,src); + std::cout< cPM; cPM(PosdefLdop,c_src); + // std::cout< ccPM; ccPM(IRLHermOpL2,cc_src); + + std::cout< +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 + +using namespace std; +using namespace Grid; +/* Params + * Grid: + * block1(4) + * block2(4) + * + * Subspace + * * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100 + * * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100 + + * Smoother: + * * Fine: Cheby(hi, lo, order) -- 60,0.5,10 + * * Coarse: Cheby(hi, lo, order) -- 12,0.1,4 + + * Lanczos: + * CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61 + */ +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template class MirsSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & SmootherMatrix; + FineOperator & SmootherOperator; + RealD tol; + RealD shift; + int maxit; + + MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : + shift(_shift),tol(_tol),maxit(_maxit), + SmootherOperator(_SmootherOperator), + SmootherMatrix(_SmootherMatrix) + {}; + + void operator() (const Field &in, Field &out) + { + ZeroGuesser Guess; + ConjugateGradient CG(tol,maxit,false); + + Field src(in.Grid()); + + ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); + SmootherOperator.AdjOp(in,src); + Guess(src,out); + CG(MdagMOp,src,out); + } +}; + +template +class MultiGridPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef CoarsenedMatrix CoarseOperator; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + + Aggregates & _Aggregates; + CoarseOperator & _CoarseOperator; + Matrix & _FineMatrix; + FineOperator & _FineOperator; + Guesser & _Guess; + FineSmoother & _Smoother; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + +#define GridLogLevel std::cout << GridLogMessage < block ({2,2,2,2}); + std::vector blockc ({2,2,2,2}); + const int nbasis= 40; + const int nbasisc= 40; + 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); + LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; + LatticeFermion result(FGrid); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + // std::string file("./ckpoint_lat.4000"); + // std::string file("./ckpoint_lat.1000"); + // NerscIO::readConfiguration(Umu,header,file); + SU::HotConfiguration(RNG4,Umu); + + 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; + LatticeFermion A(FGrid); + LatticeFermion B(FGrid); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0); + Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,400,50,50,0.0); // Slightly faster + + for(int n=0;n Level1Op; + typedef CoarsenedMatrix,nbasisc> Level2Op; + + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); + + + GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d); + std::cout << " Making 5D coarse RB grid " <,nbasisc> CoarseSubspace; + // CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< PosdefLdop(LDOp); + /* + { + int nb=nbasisc/2; + CoarseAggregates.CreateSubspaceChebyshev(CRNG,PosdefLdop,nb,15.0,0.02,1000,800,100,0.0); + for(int n=0;noSites();site++){ + subspace_g5[site](nn) = subspace[site](nn); + subspace_g5[site](nn+nb)=-subspace[site](nn+nb); + } + } + } + } + */ + typedef Level2Op::CoarseVector CoarseCoarseVector; + /* + Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix + HermitianLinearOperator L1LinOp(LDOp); + L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates); + + + std::cout< IRLHermOpL2(L2Op); + CoarseCoarseVector cc_src(CoarseCoarse5d); cc_src=1.0; + */ + /* + Chebyshev IRLChebyL2(0.001,15.0,301); + FunctionHermOp IRLOpChebyL2(IRLChebyL2,IRLHermOpL2); + PlainHermOp IRLOpL2 (IRLHermOpL2); + int cNk=24; + int cNm=36; + int cNstop=24; + ImplicitlyRestartedLanczos IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20); + + int cNconv; + std::vector eval2(cNm); + std::vector evec2(cNm,CoarseCoarse5d); + IRLL2.calc(eval2,evec2,cc_src,cNconv); + + ConjugateGradient CoarseCoarseCG(0.1,1000); + DeflatedGuesser DeflCoarseCoarseGuesser(evec2,eval2); + NormalEquations DeflCoarseCoarseCGNE(L2Op,CoarseCoarseCG,DeflCoarseCoarseGuesser); + */ + + /* + std::cout< IRLHermOp(LDOp); + // Chebyshev IRLCheby(0.001,15.0,301); + Chebyshev IRLCheby(0.03,12.0,101); + FunctionHermOp IRLOpCheby(IRLCheby,IRLHermOp); + PlainHermOp IRLOp (IRLHermOp); + int Nk=64; + int Nm=128; + int Nstop=Nk; + 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); + */ + CoarseVector c_src(Coarse5d); c_src=1.0; + // DeflatedGuesser DeflCoarseGuesser(evec,eval); + // NormalEquations DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser); + + std::cout< , NormalEquations > TwoLevelMG; + typedef MultiGridPreconditioner , NormalEquations > TwoLevelMG; + typedef MultiGridPreconditioner,nbasisc,Level1Op, DeflatedGuesser, NormalEquations > CoarseMG; + typedef MultiGridPreconditioner, LinearFunction > ThreeLevelMG; + + ChebyshevSmoother FineSmoother(0.25,60.0,12,HermIndefOp,Ddwf); + /* + // MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space + ChebyshevSmoother CoarseSmoother(0.1,15.0,3,L1LinOp,LDOp); + + // MirsSmoother CoarseCGSmoother(0.1,0.1,4,L1LinOp,LDOp); + // MirsSmoother FineCGSmoother(0.0,0.01,8,HermIndefOp,Ddwf); + + CoarseMG Level2Precon (CoarseAggregates, L2Op, + L1LinOp,LDOp, + CoarseSmoother, + DeflCoarseCoarseGuesser, + DeflCoarseCoarseCGNE); + Level2Precon.Level(2); + + // PGCR Applying this solver to solve the coarse space problem + PrecGeneralisedConjugateResidual l2PGCR(0.1, 100, L1LinOp,Level2Precon,16,16); + l2PGCR.Level(2); + + // Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space + ZeroGuesser CoarseZeroGuesser; + ThreeLevelMG ThreeLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + l2PGCR); + ThreeLevelPrecon.Level(1); + + // Apply the fine-coarse-coarsecoarse 2 deep MG preconditioner in an outer PGCR on the fine fgrid + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,1000,HermIndefOp,ThreeLevelPrecon,16,16); + l1PGCR.Level(1); + */ + std::cout< CoarseZeroGuesser; + ConjugateGradient CoarseCG(0.01,1000); + NormalEquations CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser); + TwoLevelMG TwoLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + CoarseCGNE); + TwoLevelPrecon.Level(1); + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16); + l1PGCR.Level(1); + l1PGCR(src,result); + + std::cout< pCG(1.0e-8,60000); + result=Zero(); + // pCG(HermDefOp,src,result); + + std::cout< HermOpEO(Ddwf); + // pCG(HermOpEO,src_o,result_o); + + std::cout< PM; PM(HermDefOp,src); + std::cout< cPM; cPM(PosdefLdop,c_src); + // std::cout< ccPM; ccPM(IRLHermOpL2,cc_src); + + std::cout< +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 + +using namespace std; +using namespace Grid; +/* Params + * Grid: + * block1(4) + * block2(4) + * + * Subspace + * * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100 + * * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100 + + * Smoother: + * * Fine: Cheby(hi, lo, order) -- 60,0.5,10 + * * Coarse: Cheby(hi, lo, order) -- 12,0.1,4 + + * Lanczos: + * CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61 + */ +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class SolverWrapper : public LinearFunction { +private: + CheckerBoardedSparseMatrixBase & _Matrix; + SchurRedBlackBase & _Solver; +public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations trick + ///////////////////////////////////////////////////// + SolverWrapper(CheckerBoardedSparseMatrixBase &Matrix, + SchurRedBlackBase &Solver) + : _Matrix(Matrix), _Solver(Solver) {}; + + void operator() (const Field &in, Field &out){ + + _Solver(_Matrix,in,out); // Mdag M out = Mdag in + + } +}; + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template class MirsSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & SmootherMatrix; + FineOperator & SmootherOperator; + RealD tol; + RealD shift; + int maxit; + + MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : + shift(_shift),tol(_tol),maxit(_maxit), + SmootherOperator(_SmootherOperator), + SmootherMatrix(_SmootherMatrix) + {}; + + void operator() (const Field &in, Field &out) + { + ZeroGuesser Guess; + ConjugateGradient CG(tol,maxit,false); + + Field src(in.Grid()); + + ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); + SmootherOperator.AdjOp(in,src); + Guess(src,out); + CG(MdagMOp,src,out); + } +}; + +template +class MultiGridPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef CoarsenedMatrix CoarseOperator; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + + Aggregates & _Aggregates; + CoarseOperator & _CoarseOperator; + Matrix & _FineMatrix; + FineOperator & _FineOperator; + Guesser & _Guess; + FineSmoother & _Smoother; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + +#define GridLogLevel std::cout << GridLogMessage < block ({2,2,2,2}); + //std::vector block ({2,2,2,2}); + const int nbasis= 40; + const int nbasisc= 40; + 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); + LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; + LatticeFermion result(FGrid); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + //std::string file("./ckpoint_lat.4000"); + std::string file("./ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + 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; + LatticeFermion A(FGrid); + LatticeFermion B(FGrid); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0); + Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster + + for(int n=0;n Level1Op; + typedef CoarsenedMatrix,nbasisc> Level2Op; + + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); + + + GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d); + GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d); + + Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); + + ////////////////////////////////////////////////// + // Deflate the course space. Recursive multigrid? + ////////////////////////////////////////////////// + typedef Aggregation,nbasisc> CoarseSubspace; + // CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< PosdefLdop(LDOp); + typedef Level2Op::CoarseVector CoarseCoarseVector; + CoarseVector c_src(Coarse5d); c_src=1.0; + + std::cout< , SolverWrapper > TwoLevelMG; + typedef MultiGridPreconditioner,nbasisc,Level1Op, DeflatedGuesser, NormalEquations > CoarseMG; + typedef MultiGridPreconditioner, LinearFunction > ThreeLevelMG; + + std::cout< tols({0.015}); + std::vector ords({12}); + std::vector los({0.8}); + for(int l=0;l FineSmoother(los[l],60.0,ords[o],HermIndefOp,Ddwf); + ZeroGuesser CoarseZeroGuesser; + ConjugateGradient CoarseCG(tols[t],10000); + SchurRedBlackDiagMooeeSolve CoarseRBCG(CoarseCG); + SolverWrapper CoarseSolver(LDOp,CoarseRBCG); + + TwoLevelMG TwoLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + CoarseSolver); + TwoLevelPrecon.Level(1); + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16); + l1PGCR.Level(1); + l1PGCR(src,result); + }}} + + ConjugateGradient pCG(1.0e-8,60000); + std::cout< HermOpEO(Ddwf); + pCG(HermOpEO,src_o,result_o); + + std::cout< PM; PM(HermDefOp,src); + std::cout< cPM; cPM(PosdefLdop,c_src); + // std::cout< ccPM; ccPM(IRLHermOpL2,cc_src); + + std::cout< +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 + +using namespace std; +using namespace Grid; +/* Params + * Grid: + * block1(4) + * block2(4) + * + * Subspace + * * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100 + * * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100 + + * Smoother: + * * Fine: Cheby(hi, lo, order) -- 60,0.5,10 + * * Coarse: Cheby(hi, lo, order) -- 12,0.1,4 + + * Lanczos: + * CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61 + */ +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template class MirsSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & SmootherMatrix; + FineOperator & SmootherOperator; + RealD tol; + RealD shift; + int maxit; + + MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) : + shift(_shift),tol(_tol),maxit(_maxit), + SmootherOperator(_SmootherOperator), + SmootherMatrix(_SmootherMatrix) + {}; + + void operator() (const Field &in, Field &out) + { + ZeroGuesser Guess; + ConjugateGradient CG(tol,maxit,false); + + Field src(in.Grid()); + + ShiftedMdagMLinearOperator,Field> MdagMOp(SmootherMatrix,shift); + SmootherOperator.AdjOp(in,src); + Guess(src,out); + CG(MdagMOp,src,out); + } +}; + +template +class MultiGridPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef CoarsenedMatrix CoarseOperator; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + + Aggregates & _Aggregates; + CoarseOperator & _CoarseOperator; + Matrix & _FineMatrix; + FineOperator & _FineOperator; + Guesser & _Guess; + FineSmoother & _Smoother; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + +#define GridLogLevel std::cout << GridLogMessage < block ({2,2,2,2}); + std::vector blockc ({2,2,2,2}); + const int nbasis= 40; + const int nbasisc= 40; + 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); + LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; + LatticeFermion result(FGrid); + LatticeGaugeField Umu(UGrid); + + FieldMetaData header; + // std::string file("./ckpoint_lat.4000"); + std::string file("./ckpoint_lat.1000"); + NerscIO::readConfiguration(Umu,header,file); + + 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; + LatticeFermion A(FGrid); + LatticeFermion B(FGrid); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0); + // Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0); + Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster + + for(int n=0;n Level1Op; + typedef CoarsenedMatrix,nbasisc> Level2Op; + + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); + + + GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d); + std::cout << " Making 5D coarse RB grid " <,nbasisc> CoarseSubspace; + // CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< PosdefLdop(LDOp); + /* + { + int nb=nbasisc/2; + CoarseAggregates.CreateSubspaceChebyshev(CRNG,PosdefLdop,nb,15.0,0.02,1000,800,100,0.0); + for(int n=0;noSites();site++){ + subspace_g5[site](nn) = subspace[site](nn); + subspace_g5[site](nn+nb)=-subspace[site](nn+nb); + } + } + } + } + */ + typedef Level2Op::CoarseVector CoarseCoarseVector; + /* + Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix + HermitianLinearOperator L1LinOp(LDOp); + L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates); + + + std::cout< IRLHermOpL2(L2Op); + CoarseCoarseVector cc_src(CoarseCoarse5d); cc_src=1.0; + */ + /* + Chebyshev IRLChebyL2(0.001,15.0,301); + FunctionHermOp IRLOpChebyL2(IRLChebyL2,IRLHermOpL2); + PlainHermOp IRLOpL2 (IRLHermOpL2); + int cNk=24; + int cNm=36; + int cNstop=24; + ImplicitlyRestartedLanczos IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20); + + int cNconv; + std::vector eval2(cNm); + std::vector evec2(cNm,CoarseCoarse5d); + IRLL2.calc(eval2,evec2,cc_src,cNconv); + + ConjugateGradient CoarseCoarseCG(0.1,1000); + DeflatedGuesser DeflCoarseCoarseGuesser(evec2,eval2); + NormalEquations DeflCoarseCoarseCGNE(L2Op,CoarseCoarseCG,DeflCoarseCoarseGuesser); + */ + + /* + std::cout< IRLHermOp(LDOp); + // Chebyshev IRLCheby(0.001,15.0,301); + Chebyshev IRLCheby(0.03,12.0,101); + FunctionHermOp IRLOpCheby(IRLCheby,IRLHermOp); + PlainHermOp IRLOp (IRLHermOp); + int Nk=64; + int Nm=128; + int Nstop=Nk; + 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); + */ + CoarseVector c_src(Coarse5d); c_src=1.0; + // DeflatedGuesser DeflCoarseGuesser(evec,eval); + // NormalEquations DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser); + + std::cout< , NormalEquations > TwoLevelMG; + typedef MultiGridPreconditioner , NormalEquations > TwoLevelMG; + typedef MultiGridPreconditioner,nbasisc,Level1Op, DeflatedGuesser, NormalEquations > CoarseMG; + typedef MultiGridPreconditioner, LinearFunction > ThreeLevelMG; + + ChebyshevSmoother FineSmoother(0.25,60.0,12,HermIndefOp,Ddwf); + /* + // MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space + ChebyshevSmoother CoarseSmoother(0.1,15.0,3,L1LinOp,LDOp); + + // MirsSmoother CoarseCGSmoother(0.1,0.1,4,L1LinOp,LDOp); + // MirsSmoother FineCGSmoother(0.0,0.01,8,HermIndefOp,Ddwf); + + CoarseMG Level2Precon (CoarseAggregates, L2Op, + L1LinOp,LDOp, + CoarseSmoother, + DeflCoarseCoarseGuesser, + DeflCoarseCoarseCGNE); + Level2Precon.Level(2); + + // PGCR Applying this solver to solve the coarse space problem + PrecGeneralisedConjugateResidual l2PGCR(0.1, 100, L1LinOp,Level2Precon,16,16); + l2PGCR.Level(2); + + // Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space + ZeroGuesser CoarseZeroGuesser; + ThreeLevelMG ThreeLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + l2PGCR); + ThreeLevelPrecon.Level(1); + + // Apply the fine-coarse-coarsecoarse 2 deep MG preconditioner in an outer PGCR on the fine fgrid + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,1000,HermIndefOp,ThreeLevelPrecon,16,16); + l1PGCR.Level(1); + */ + std::cout< CoarseZeroGuesser; + ConjugateGradient CoarseCG(0.01,1000); + NormalEquations CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser); + TwoLevelMG TwoLevelPrecon(Aggregates, LDOp, + HermIndefOp,Ddwf, + FineSmoother, + CoarseZeroGuesser, + CoarseCGNE); + TwoLevelPrecon.Level(1); + PrecGeneralisedConjugateResidual l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16); + l1PGCR.Level(1); + l1PGCR(src,result); + + std::cout< pCG(1.0e-8,60000); + result=Zero(); + // pCG(HermDefOp,src,result); + + std::cout< HermOpEO(Ddwf); + pCG(HermOpEO,src_o,result_o); + + std::cout< PM; PM(HermDefOp,src); + std::cout< cPM; cPM(PosdefLdop,c_src); + // std::cout< ccPM; ccPM(IRLHermOpL2,cc_src); + + std::cout< + + 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; + +// TODO +// +// Coarse Grid axpby_ssp_pminus // Inherit from spProj5pm +// Coarse Grid axpby_ssp_pplus + +template +class CayleyBase : public SparseMatrixBase +{ +public: + int Ls; + // protected: + RealD mass; + RealD M5; + // Save arguments to SetCoefficientsInternal + Vector _gamma; + RealD _zolo_hi; + RealD _b; + RealD _c; + + // Cayley form Moebius (tanh and zolotarev) + Vector omega; + Vector bs; // S dependent coeffs + Vector cs; + Vector as; + // For preconditioning Cayley form + Vector bee; + Vector cee; + Vector aee; + Vector beo; + Vector ceo; + Vector aeo; + // LDU factorisation of the eeoo matrix + Vector lee; + Vector leem; + Vector uee; + Vector ueem; + Vector dee; +public: + CayleyBase(RealD _M5, RealD _mass, int _Ls, RealD b_, RealD c_) : + M5(_M5), + mass(_mass), + Ls(_Ls), + _b(b_), + _c(c_) + { + RealD eps = 1.0; + Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham + this->SetCoefficientsTanh(zdata,1.0,0.0); + Approx::zolotarev_free(zdata); + } + ///////////////////////////////////////////////////////// + // Replicates functionality + // Use a common base class approach + ///////////////////////////////////////////////////////// + // Tanh + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c) + { + Vector gamma(this->Ls); + for(int s=0;sLs;s++) gamma[s] = zdata->gamma[s]; + SetCoefficientsInternal(1.0,gamma,b,c); + } + //Zolo + void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c) + { + Vector gamma(this->Ls); + for(int s=0;sLs;s++) gamma[s] = zdata->gamma[s]; + SetCoefficientsInternal(zolo_hi,gamma,b,c); + } + //Zolo + void SetCoefficientsInternal(RealD zolo_hi,Vector & gamma,RealD b,RealD c) + { + int Ls=this->Ls; + + /////////////////////////////////////////////////////////// + // The Cayley coeffs (unprec) + /////////////////////////////////////////////////////////// + assert(gamma.size()==Ls); + + omega.resize(Ls); + bs.resize(Ls); + cs.resize(Ls); + as.resize(Ls); + + double bpc = b+c; + double bmc = b-c; + _b = b; + _c = c; + _gamma = gamma; // Save the parameters so we can change mass later. + _zolo_hi= zolo_hi; + for(int i=0; i < Ls; i++){ + as[i] = 1.0; + omega[i] = _gamma[i]*_zolo_hi; //NB reciprocal relative to Chroma NEF code + assert(omega[i]!=Coeff_t(0.0)); + bs[i] = 0.5*(bpc/omega[i] + bmc); + cs[i] = 0.5*(bpc/omega[i] - bmc); + } + + //////////////////////////////////////////////////////// + // Constants for the preconditioned matrix Cayley form + //////////////////////////////////////////////////////// + bee.resize(Ls); + cee.resize(Ls); + beo.resize(Ls); + ceo.resize(Ls); + + for(int i=0;iM5) +1.0); + assert(bee[i]!=Coeff_t(0.0)); + cee[i]=as[i]*(1.0-cs[i]*(4.0-this->M5)); + beo[i]=as[i]*bs[i]; + ceo[i]=-as[i]*cs[i]; + } + aee.resize(Ls); + aeo.resize(Ls); + for(int i=0;i &out){assert(0);}; + virtual void DW (const Field &psi, Field &chi)=0; + virtual void DWDag (const Field &psi, Field &chi)=0; + + void M (const Field &psi, Field &chi) + { + Field Din(psi.Grid()); + Meooe5D(psi,Din); + DW(Din,chi); + axpby(chi,1.0,1.0,chi,psi); + M5D(psi,chi); + } + void Mdag (const Field &psi, Field &chi) + { + Field Din(psi.Grid()); + DWDag(psi,Din); + MeooeDag5D(Din,chi); + M5Ddag(psi,chi); + axpby (chi,1.0,1.0,chi,psi); + } + ///////////////////////////////// + // P and Pdag - might be needed + ///////////////////////////////// + void P(const Field &psi, Field &chi) + { + int Ls= this->Ls; + chi=Zero(); + for(int s=0;sLs; + chi=Zero(); + for(int s=0;sLs; + Vector diag (Ls,1.0); + Vector upper(Ls,-1.0); upper[Ls-1]=mass; + Vector lower(Ls,-1.0); lower[0] =mass; + M5D(psi,chi,chi,lower,diag,upper); + } + void M5Ddag (const Field &psi, Field &chi) + { + int Ls=this->Ls; + Vector diag(Ls,1.0); + Vector upper(Ls,-1.0); + Vector lower(Ls,-1.0); + upper[Ls-1]=-mass*upper[Ls-1]; + lower[0] =-mass*lower[0]; + M5Ddag(psi,chi,chi,lower,diag,upper); + } + void Meooe5D (const Field &psi, Field &Din) + { + int Ls=this->Ls; + Vector diag = bs; + Vector upper= cs; + Vector lower= cs; + upper[Ls-1]=-mass*upper[Ls-1]; + lower[0] =-mass*lower[0]; + M5D(psi,psi,Din,lower,diag,upper); + } + void MeooeDag5D (const Field &psi, Field &Din) + { + int Ls=this->Ls; + Vector diag =bs; + Vector upper=cs; + Vector lower=cs; + + for (int s=0;s &lower, + Vector &diag, + Vector &upper) + { + chi_i.Checkerboard()=psi_i.Checkerboard(); + GridBase *grid=psi_i.Grid(); + autoView(psi , psi_i,AcceleratorRead); + autoView(phi , phi_i,AcceleratorRead); + autoView(chi , chi_i,AcceleratorWrite); + assert(phi.Checkerboard() == psi.Checkerboard()); + + auto pdiag = &diag[0]; + auto pupper = &upper[0]; + auto plower = &lower[0]; + + int Ls =this->Ls; + + // 10 = 3 complex mult + 2 complex add + // Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting) + uint64_t nloop = grid->oSites()/Ls; + + const int Nsimd = Field::vector_type::Nsimd(); + accelerator_for(sss,nloop,Nsimd,{ + uint64_t ss= sss*Ls; + typedef decltype(coalescedRead(psi[0])) spinor; + spinor tmp1, tmp2; + for(int s=0;s &lower, + Vector &diag, + Vector &upper) + { + chi_i.Checkerboard()=psi_i.Checkerboard(); + GridBase *grid=psi_i.Grid(); + autoView(psi , psi_i,AcceleratorRead); + autoView(phi , phi_i,AcceleratorRead); + autoView(chi , chi_i,AcceleratorWrite); + assert(phi.Checkerboard() == psi.Checkerboard()); + + auto pdiag = &diag[0]; + auto pupper = &upper[0]; + auto plower = &lower[0]; + + int Ls=this->Ls; + + uint64_t nloop = grid->oSites()/Ls; + const int Nsimd = Field::vector_type::Nsimd(); + accelerator_for(sss,nloop,Nsimd,{ + uint64_t ss=sss*Ls; + typedef decltype(coalescedRead(psi[0])) spinor; + spinor tmp1,tmp2; + for(int s=0;s +class CoarseCayleyFermion : public CayleyBase< Lattice > , ComplexD > +{ +public: + typedef iVector siteVector; + typedef Lattice CoarseComplexField; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; + typedef iMatrix Cobj; + typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field + typedef Lattice FineField; + + // Similar to the CoarseOperator but add 5D support. + Geometry geom; + GridBase *Coarse5D; + GridBase *Coarse4D; + CartesianStencil Stencil; + CoarsenedMatrix &Dw; + + GridBase * Grid(void) { return Coarse5D; }; // this is all the linalg routines need to know + + CoarseCayleyFermion(GridCartesian &CoarseGrid4, + GridCartesian &CoarseGrid5, + CoarsenedMatrix &_Dw, + RealD M5, RealD mass, int Ls, RealD b, RealD c) : + CayleyBase(M5,mass,Ls,b,c), + Coarse4D(&CoarseGrid4), + Coarse5D(&CoarseGrid5), + Dw(_Dw), + geom(CoarseGrid5._ndimension), + Stencil( &CoarseGrid5,geom.npoint,Even,geom.directions,geom.displacements,0) + { + }; + +public: + void Project( CoarseVector &C ) + { + const int Nsimd = CComplex::Nsimd(); + autoView(Cv,C, AcceleratorWrite); + int Ls = this->Ls; + for(int s=0;soSites(), Nsimd, { + int sF= sU*Ls+s; + auto tmp = coalescedRead(Cv[sF]); + coalescedWrite(Cv[sF],tmp); + }); + } + } + //////////////////////////////////////////////// + // This is specific to Coarse Grid Cayley + //////////////////////////////////////////////// + virtual void Mdiag (const CoarseVector &in, CoarseVector &out) + { + std::vector allout(9,in.Grid()); + this->MdirAll(in,allout); + out = allout[8]; + } + virtual void Mdir (const CoarseVector &in, CoarseVector &out,int dir, int disp) + { + assert(0); + } + virtual void MdirAll (const CoarseVector &in, std::vector &out) + { + conformable(Coarse5D,in.Grid()); + + SimpleCompressor compressor; + + Stencil.HaloExchange(in,compressor); + typedef LatticeView Aview; + + const int Nsimd = CComplex::Nsimd(); + + // Ls loop for2D + int Ls=this->Ls; + + siteVector *CBp=Stencil.CommBuf(); + + int ptype; + int nb2=nbasis/2; + + autoView(in_v , in, AcceleratorRead); + autoView(st, Stencil, AcceleratorRead); + for(int point=0;pointoSites(), b, nbasis, Nsimd, { + + typedef decltype(coalescedRead(in_v[0])) calcVector; + typedef decltype(coalescedRead(in_v[0](0))) calcComplex; + int sU = sF/Ls; + int s = sF%Ls; + + calcComplex res = Zero(); + calcVector nbr; + int ptype; + + StencilEntry *SE=st.GetEntry(ptype,point,sF); + + if(SE->_is_local) { + nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); + } else { + nbr = coalescedRead(CBp[SE->_offset]); + } + acceleratorSynchronise(); + + for(int bb=0;bb compressor; + + Stencil.HaloExchange(in,compressor); + typedef LatticeView Aview; + + const int Nsimd = CComplex::Nsimd(); + + // Ls loop for2D + int Ls=this->Ls; + + Vector AcceleratorViewContainer; + for(int p=0;poSites(), b, nbasis, Nsimd, { + + typedef decltype(coalescedRead(in_v[0])) calcVector; + typedef decltype(coalescedRead(in_v[0](0))) calcComplex; + int sU = sF/Ls; + int s = sF%Ls; + + calcComplex res = Zero(); + + { + calcVector nbr; + int ptype; + + for(int point=0;point_is_local) { + nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); + } else { + nbr = coalescedRead(CBp[SE->_offset]); + } + acceleratorSynchronise(); + + for(int bb=0;bb Aggregates; + + void PromoteFromSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(C4,C,s,0); + _Aggregates.PromoteFromSubspace(C4,F4); + InsertSlice(F4,F,s,0); + } + } + void ProjectToSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(F4,F,s,0); + _Aggregates.ProjectToSubspace (C4,F4); + InsertSlice(C4,C,s,0); + } + Project(C); + } + template + void Test(Aggregates &_Aggregates,GridBase *FineGrid, Ddwf &_Ddwf) + { + typedef Lattice FineField; + CoarseVector Cin(Coarse5D); + CoarseVector Cout(Coarse5D); + CoarseVector CFout(Coarse5D); + + FineField Fin(FineGrid); + FineField Fout(FineGrid); + + + std::vector seeds({1,2,3,4,5}); + GridParallelRNG RNG(Coarse5D); RNG.SeedFixedIntegers(seeds); + + gaussian(RNG,Cin); + PromoteFromSubspace(_Aggregates,Cin,Fin); + ProjectToSubspace(_Aggregates,Cin,Fin); + + std::cout << GridLogMessage<< "************ "<M(Cin,Cout); + this->Project(Cout); + std::cout << GridLogMessage<< " Cout "<Mdag(Cin,Cout); + this->Project(Cout); + std::cout << GridLogMessage<< " Cout "< Directions(void) { return geom.directions;}; + virtual std::vector Displacements(void){ return geom.displacements;}; +}; + + +template class SolverWrapper : public LinearFunction { +private: + LinearOperatorBase & _Matrix; + OperatorFunction & _Solver; + LinearFunction & _Guess; +public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations trick + ///////////////////////////////////////////////////// + SolverWrapper(LinearOperatorBase &Matrix, + OperatorFunction &Solver, + LinearFunction &Guess) + : _Matrix(Matrix), _Solver(Solver), _Guess(Guess) {}; + + void operator() (const Field &in, Field &out){ + + _Guess(in,out); + _Solver(_Matrix,in,out); // Mdag M out = Mdag in + + } +}; + +// Must use a non-hermitian solver +template +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){}; + + virtual std::vector Directions(void) { return _Mat.Directions();}; + virtual std::vector Displacements(void){ return _Mat.Displacements();}; + + 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){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(tmp,out); + _Mat.Mdag(in,tmp); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + assert(0); + } + void HermOp(const Field &in, Field &out){ + assert(0); + } +}; + +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template +class MGPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef CoarseCayleyFermion CoarseOperator; + // typedef SparseMatrixBase CoarseOperator; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + auto CoarseGrid = _CoarseOperator.Grid(); + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + std::cout< +class HDCRPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + //typedef CoarseCayleyFermion CoarseOperator; + typedef SparseMatrixBase CoarseOperator; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + + HDCRPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + auto CoarseGrid = _CoarseOperator.Grid(); + CoarseVector Csrc(CoarseGrid); + CoarseVector g5Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + std::cout< block ({2,2,2,2}); // 4,2,2,2 gets worse + std::vector blockc ({1,1,1,1}); + const int nbasis= 24; + const int nbasisc= 32; // decrease, not improvement + + auto clatt = GridDefaultLatt(); + for(int d=0;d seeds({1,2,3,4}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(seeds); + + LatticeGaugeField Umu(UGrid); +#if 0 + SU3::TepidConfiguration(RNG4,Umu); + RealD M5=1.0; +#else + std::string file("./ckpoint_lat.1000"); + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,file); + RealD M5=1.8; +#endif + + std::cout< Subspace; + typedef CoarsenedMatrix CoarseOperator; + typedef CoarseOperator::CoarseVector CoarseVector; + typedef CoarseOperator::siteVector siteVector; + + std::cout< MdagM_Dw(Dw_null); + + std::cout< WilsonCG(1.0e-10,40000); + LatticeFermion w_src(UGrid); w_src=1.0; + LatticeFermion w_res(UGrid); + WilsonCG(MdagM_Dw,w_src,w_res); + exit(0); + */ + std::cout< Level1Op4; + typedef CoarseCayleyFermion Level1Op5; + Level1Op4 c_Dw (*Coarse4d,0); + NonHermitianLinearOperator LinOpDw(Dw); + c_Dw.CoarsenOperator(UGrid,LinOpDw,Aggregates4D); // contains the M5 from Dw(-M5) + // c_Dw.Test(Aggregates4D,UGrid,LinOpDw); + + std::cout< MdagM_cDwf(c_Dwf); + + std::cout<,nbasisc> Level2Op; + typedef Aggregation,nbasisc> CoarseSubspace; + CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< L1Hdwf(c_Dwf); + GridRedBlackCartesian * CoarseCoarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,CoarseCoarse4d); + Level2Op cc_Dwf (*CoarseCoarse5d,*CoarseCoarse5dRB,1); // say it is hermitian + cc_Dwf.CoarsenOperator(Coarse5d,L1Hdwf,CoarseAggregates); + // cc_Dwf.Test(CoarseAggregates,Coarse5d,L1Hdwf); + + typedef Level2Op::CoarseVector CoarseCoarseVector; + + std::cout< CoarseCG(tol,MaxIt); + ConjugateGradient FineCG(tol,MaxIt); + + NonHermitianLinearOperator FineM(Ddwf); + MdagMLinearOperator FineMdagM(Ddwf); // M^\dag M + + NonHermitianLinearOperator CoarseM(c_Dwf); + MdagMLinearOperator CoarseMdagM(c_Dwf); + + NonHermitianLinearOperator CoarseCoarseM(cc_Dwf); + MdagMLinearOperator CoarseCoarseMdagM(cc_Dwf); + + + std::cout< PM; PM(MdagM_Dw,w_src); + std::cout< cPM; cPM(CoarseMdagM,c_src); + + cc_src=1.0; + PowerMethod ccPM; ccPM(CoarseCoarseMdagM,cc_src); + + std::cout< IRLHermOpL2(cc_Dwf); + Chebyshev IRLChebyL2(IRL_lo,IRL_hi,IRL_ord); + FunctionHermOp IRLOpChebyL2(IRLChebyL2,IRLHermOpL2); + PlainHermOp IRLOpL2 (IRLHermOpL2); + ImplicitlyRestartedLanczos IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20); + + int cNconv; + cNm=0; + std::vector eval2(cNm); + std::vector evec2(cNm,CoarseCoarse5d); + cc_src=1.0; + // IRLL2.calc(eval2,evec2,cc_src,cNconv); + + ConjugateGradient CoarseCoarseCG(0.02,10000); + DeflatedGuesser DeflCoarseCoarseGuesser(evec2,eval2); + NormalEquations DeflCoarseCoarseCGNE(cc_Dwf,CoarseCoarseCG,DeflCoarseCoarseGuesser); + + ZeroGuesser CoarseZeroGuesser; + ZeroGuesser CoarseCoarseZeroGuesser; + + std::cout< CoarseCoarseCGNE(cc_Dwf,CoarseCoarseCG,CoarseCoarseZeroGuesser); + { +typedef HDCRPreconditioner,nbasisc,NormalEquations > CoarseMG; + typedef MGPreconditioner > ThreeLevelMG; + + // MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space + ChebyshevSmoother CoarseSmoother1(0.5,22.0,12,CoarseM,c_Dwf); // 37s, 26 iter + ChebyshevSmoother CoarseSmoother2(0.5,22.0,12,CoarseM,c_Dwf); + + // ChebyshevSmoother CoarseSmoother1(0.5,22.0,7,CoarseM,c_Dwf); // 38s, 26 iter + // ChebyshevSmoother CoarseSmoother2(0.5,22.0,7,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.4,22.0,7,CoarseM,c_Dwf); // 41s, 27 iter + // ChebyshevSmoother CoarseSmoother2(0.4,22.0,7,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.6,22.0,6,CoarseM,c_Dwf); // 26 iter + // ChebyshevSmoother CoarseSmoother2(0.6,22.0,6,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.5,22.0,5,CoarseM,c_Dwf); // 33 iter, 55s + // ChebyshevSmoother CoarseSmoother2(0.5,22.0,5,CoarseM,c_Dwf); + + + CoarseMG Level2Precon (CoarseAggregates, + CoarseM, + CoarseSmoother1, + CoarseSmoother2, + cc_Dwf, + DeflCoarseCoarseCGNE); + Level2Precon.Level(2); + + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.5, 100, CoarseM,Level2Precon,16,16); // 26 iter, 37s + // PGCR Applying this solver to solve the coarse space problem + // COULD BE FIXED??? + PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.0, 1, CoarseM,Level2Precon,2,2); + + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(1.0, 100, CoarseM,Level2Precon,16,16); // 35 iter, 45s + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.6, 100, CoarseM,Level2Precon,16,16); // 26,38 (diifferene is measurement noise) + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.2, 100, CoarseM,Level2Precon,16,16); // 26 iter, 47s + L2PGCR.Level(2); + + // Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space + + // ChebyshevSmoother FineSmoother1(0.5,60.0,14,FineM,Ddwf); // 26 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,14,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 25 iter, 38s + // ChebyshevSmoother FineSmoother2(0.5,60.0,16,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 23 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,20,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,10,FineM,Ddwf);24 iter, 44s + // ChebyshevSmoother FineSmoother2(0.5,60.0,24,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // odd convergence tail at 10^-9 ish + // ChebyshevSmoother FineSmoother2(0.1,60.0,24,FineM,Ddwf); // 33 iter, waas O(10-9 by 26) + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 25 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,18,FineM,Ddwf); // + + ChebyshevSmoother FineSmoother1(0.5,60.0,16,FineM,Ddwf); + ChebyshevSmoother FineSmoother2(0.5,60.0,16,FineM,Ddwf); // + + // ChebyshevSmoother FineSmoother1(0.5,60.0,11,FineM,Ddwf); // 33 iter, 49s + // ChebyshevSmoother FineSmoother2(0.5,60.0,11,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 26 iter, 37s + // ChebyshevSmoother FineSmoother2(0.5,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.4,60.0,12,FineM,Ddwf); // iter 26 no change in final residual + // ChebyshevSmoother FineSmoother2(0.4,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.3,60.0,12,FineM,Ddwf); // 27 iter 39s. + // ChebyshevSmoother FineSmoother2(0.3,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.3,60.0,13,FineM,Ddwf); // 26 iter, but slower + // ChebyshevSmoother FineSmoother2(0.3,60.0,13,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(1.0,60.0,12,FineM,Ddwf); // 34 iter, slower + // ChebyshevSmoother FineSmoother2(1.0,60.0,12,FineM,Ddwf); + + ThreeLevelMG ThreeLevelPrecon(Aggregates4D, + FineM, + FineSmoother1, + FineSmoother2, + c_Dwf, + L2PGCR); + ThreeLevelPrecon.Level(1); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,FineM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + + f_res=Zero(); + L1PGCR(f_src,f_res); + } + + std::cout< + + 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; + +// TODO +// +// Coarse Grid axpby_ssp_pminus // Inherit from spProj5pm +// Coarse Grid axpby_ssp_pplus + +template +class CayleyBase : public SparseMatrixBase +{ +public: + int Ls; + // protected: + RealD mass; + RealD M5; + // Save arguments to SetCoefficientsInternal + Vector _gamma; + RealD _zolo_hi; + RealD _b; + RealD _c; + + // Cayley form Moebius (tanh and zolotarev) + Vector omega; + Vector bs; // S dependent coeffs + Vector cs; + Vector as; + // For preconditioning Cayley form + Vector bee; + Vector cee; + Vector aee; + Vector beo; + Vector ceo; + Vector aeo; + // LDU factorisation of the eeoo matrix + Vector lee; + Vector leem; + Vector uee; + Vector ueem; + Vector dee; +public: + CayleyBase(RealD _M5, RealD _mass, int _Ls, RealD b_, RealD c_) : + M5(_M5), + mass(_mass), + Ls(_Ls), + _b(b_), + _c(c_) + { + RealD eps = 1.0; + Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham + this->SetCoefficientsTanh(zdata,1.0,0.0); + Approx::zolotarev_free(zdata); + } + ///////////////////////////////////////////////////////// + // Replicates functionality + // Use a common base class approach + ///////////////////////////////////////////////////////// + // Tanh + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c) + { + Vector gamma(this->Ls); + for(int s=0;sLs;s++) gamma[s] = zdata->gamma[s]; + SetCoefficientsInternal(1.0,gamma,b,c); + } + //Zolo + void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c) + { + Vector gamma(this->Ls); + for(int s=0;sLs;s++) gamma[s] = zdata->gamma[s]; + SetCoefficientsInternal(zolo_hi,gamma,b,c); + } + //Zolo + void SetCoefficientsInternal(RealD zolo_hi,Vector & gamma,RealD b,RealD c) + { + int Ls=this->Ls; + + /////////////////////////////////////////////////////////// + // The Cayley coeffs (unprec) + /////////////////////////////////////////////////////////// + assert(gamma.size()==Ls); + + omega.resize(Ls); + bs.resize(Ls); + cs.resize(Ls); + as.resize(Ls); + + double bpc = b+c; + double bmc = b-c; + _b = b; + _c = c; + _gamma = gamma; // Save the parameters so we can change mass later. + _zolo_hi= zolo_hi; + for(int i=0; i < Ls; i++){ + as[i] = 1.0; + omega[i] = _gamma[i]*_zolo_hi; //NB reciprocal relative to Chroma NEF code + assert(omega[i]!=Coeff_t(0.0)); + bs[i] = 0.5*(bpc/omega[i] + bmc); + cs[i] = 0.5*(bpc/omega[i] - bmc); + } + + //////////////////////////////////////////////////////// + // Constants for the preconditioned matrix Cayley form + //////////////////////////////////////////////////////// + bee.resize(Ls); + cee.resize(Ls); + beo.resize(Ls); + ceo.resize(Ls); + + for(int i=0;iM5) +1.0); + assert(bee[i]!=Coeff_t(0.0)); + cee[i]=as[i]*(1.0-cs[i]*(4.0-this->M5)); + beo[i]=as[i]*bs[i]; + ceo[i]=-as[i]*cs[i]; + } + aee.resize(Ls); + aeo.resize(Ls); + for(int i=0;i &out){assert(0);}; + virtual void DW (const Field &psi, Field &chi)=0; + virtual void DWDag (const Field &psi, Field &chi)=0; + + void M (const Field &psi, Field &chi) + { + Field Din(psi.Grid()); + Meooe5D(psi,Din); + DW(Din,chi); + axpby(chi,1.0,1.0,chi,psi); + M5D(psi,chi); + } + void Mdag (const Field &psi, Field &chi) + { + Field Din(psi.Grid()); + DWDag(psi,Din); + MeooeDag5D(Din,chi); + M5Ddag(psi,chi); + axpby (chi,1.0,1.0,chi,psi); + } + ///////////////////////////////// + // P and Pdag - might be needed + ///////////////////////////////// + void P(const Field &psi, Field &chi) + { + int Ls= this->Ls; + chi=Zero(); + for(int s=0;sLs; + chi=Zero(); + for(int s=0;sLs; + Vector diag (Ls,1.0); + Vector upper(Ls,-1.0); upper[Ls-1]=mass; + Vector lower(Ls,-1.0); lower[0] =mass; + M5D(psi,chi,chi,lower,diag,upper); + } + void M5Ddag (const Field &psi, Field &chi) + { + int Ls=this->Ls; + Vector diag(Ls,1.0); + Vector upper(Ls,-1.0); + Vector lower(Ls,-1.0); + upper[Ls-1]=-mass*upper[Ls-1]; + lower[0] =-mass*lower[0]; + M5Ddag(psi,chi,chi,lower,diag,upper); + } + void Meooe5D (const Field &psi, Field &Din) + { + int Ls=this->Ls; + Vector diag = bs; + Vector upper= cs; + Vector lower= cs; + upper[Ls-1]=-mass*upper[Ls-1]; + lower[0] =-mass*lower[0]; + M5D(psi,psi,Din,lower,diag,upper); + } + void MeooeDag5D (const Field &psi, Field &Din) + { + int Ls=this->Ls; + Vector diag =bs; + Vector upper=cs; + Vector lower=cs; + + for (int s=0;s &lower, + Vector &diag, + Vector &upper) + { + chi_i.Checkerboard()=psi_i.Checkerboard(); + GridBase *grid=psi_i.Grid(); + autoView(psi , psi_i,AcceleratorRead); + autoView(phi , phi_i,AcceleratorRead); + autoView(chi , chi_i,AcceleratorWrite); + assert(phi.Checkerboard() == psi.Checkerboard()); + + auto pdiag = &diag[0]; + auto pupper = &upper[0]; + auto plower = &lower[0]; + + int Ls =this->Ls; + + // 10 = 3 complex mult + 2 complex add + // Flops = 10.0*(Nc*Ns) *Ls*vol (/2 for red black counting) + uint64_t nloop = grid->oSites()/Ls; + + const int Nsimd = Field::vector_type::Nsimd(); + accelerator_for(sss,nloop,Nsimd,{ + uint64_t ss= sss*Ls; + typedef decltype(coalescedRead(psi[0])) spinor; + spinor tmp1, tmp2; + for(int s=0;s &lower, + Vector &diag, + Vector &upper) + { + chi_i.Checkerboard()=psi_i.Checkerboard(); + GridBase *grid=psi_i.Grid(); + autoView(psi , psi_i,AcceleratorRead); + autoView(phi , phi_i,AcceleratorRead); + autoView(chi , chi_i,AcceleratorWrite); + assert(phi.Checkerboard() == psi.Checkerboard()); + + auto pdiag = &diag[0]; + auto pupper = &upper[0]; + auto plower = &lower[0]; + + int Ls=this->Ls; + + uint64_t nloop = grid->oSites()/Ls; + const int Nsimd = Field::vector_type::Nsimd(); + accelerator_for(sss,nloop,Nsimd,{ + uint64_t ss=sss*Ls; + typedef decltype(coalescedRead(psi[0])) spinor; + spinor tmp1,tmp2; + for(int s=0;s +class CoarseCayleyFermion : public CayleyBase< Lattice > , ComplexD > +{ +public: + typedef iVector siteVector; + typedef Lattice CoarseComplexField; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; + typedef iMatrix Cobj; + typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field + typedef Lattice FineField; + + // Similar to the CoarseOperator but add 5D support. + Geometry geom; + GridBase *Coarse5D; + GridBase *Coarse4D; + CartesianStencil Stencil; + CoarsenedMatrix &Dw; + + GridBase * Grid(void) { return Coarse5D; }; // this is all the linalg routines need to know + + CoarseCayleyFermion(GridCartesian &CoarseGrid4, + GridCartesian &CoarseGrid5, + CoarsenedMatrix &_Dw, + RealD M5, RealD mass, int Ls, RealD b, RealD c) : + CayleyBase(M5,mass,Ls,b,c), + Coarse4D(&CoarseGrid4), + Coarse5D(&CoarseGrid5), + Dw(_Dw), + geom(CoarseGrid5._ndimension), + Stencil( &CoarseGrid5,geom.npoint,Even,geom.directions,geom.displacements,0) + { + }; + +public: + void Project( CoarseVector &C ) + { + const int Nsimd = CComplex::Nsimd(); + autoView(Cv,C, AcceleratorWrite); + int Ls = this->Ls; + for(int s=0;soSites(), Nsimd, { + int sF= sU*Ls+s; + auto tmp = coalescedRead(Cv[sF]); + coalescedWrite(Cv[sF],tmp); + }); + } + } + //////////////////////////////////////////////// + // This is specific to Coarse Grid Cayley + //////////////////////////////////////////////// + virtual void Mdiag (const CoarseVector &in, CoarseVector &out) + { + std::vector allout(9,in.Grid()); + this->MdirAll(in,allout); + out = allout[8]; + } + virtual void Mdir (const CoarseVector &in, CoarseVector &out,int dir, int disp) + { + assert(0); + } + virtual void MdirAll (const CoarseVector &in, std::vector &out) + { + conformable(Coarse5D,in.Grid()); + + SimpleCompressor compressor; + + Stencil.HaloExchange(in,compressor); + typedef LatticeView Aview; + + const int Nsimd = CComplex::Nsimd(); + + // Ls loop for2D + int Ls=this->Ls; + + siteVector *CBp=Stencil.CommBuf(); + + int ptype; + int nb2=nbasis/2; + + autoView(in_v , in, AcceleratorRead); + autoView(st, Stencil, AcceleratorRead); + for(int point=0;pointoSites(), b, nbasis, Nsimd, { + + typedef decltype(coalescedRead(in_v[0])) calcVector; + typedef decltype(coalescedRead(in_v[0](0))) calcComplex; + int sU = sF/Ls; + int s = sF%Ls; + + calcComplex res = Zero(); + calcVector nbr; + int ptype; + + StencilEntry *SE=st.GetEntry(ptype,point,sF); + + if(SE->_is_local) { + nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); + } else { + nbr = coalescedRead(CBp[SE->_offset]); + } + acceleratorSynchronise(); + + for(int bb=0;bb compressor; + + Stencil.HaloExchange(in,compressor); + typedef LatticeView Aview; + + const int Nsimd = CComplex::Nsimd(); + + // Ls loop for2D + int Ls=this->Ls; + + Vector AcceleratorViewContainer; + for(int p=0;poSites(), b, nbasis, Nsimd, { + + typedef decltype(coalescedRead(in_v[0])) calcVector; + typedef decltype(coalescedRead(in_v[0](0))) calcComplex; + int sU = sF/Ls; + int s = sF%Ls; + + calcComplex res = Zero(); + + { + calcVector nbr; + int ptype; + + for(int point=0;point_is_local) { + nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); + } else { + nbr = coalescedRead(CBp[SE->_offset]); + } + acceleratorSynchronise(); + + for(int bb=0;bb Aggregates; + + void PromoteFromSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(C4,C,s,0); + _Aggregates.PromoteFromSubspace(C4,F4); + InsertSlice(F4,F,s,0); + } + } + void ProjectToSubspace(Aggregates &_Aggregates,CoarseVector &C,FineField &F) + { + auto FineGrid4 = _Aggregates.FineGrid; + FineField F4(FineGrid4); + CoarseVector C4(Coarse4D); + for(int s=0;sLs;s++){ + ExtractSlice(F4,F,s,0); + _Aggregates.ProjectToSubspace (C4,F4); + InsertSlice(C4,C,s,0); + } + Project(C); + } + template + void Test(Aggregates &_Aggregates,GridBase *FineGrid, Ddwf &_Ddwf) + { + typedef Lattice FineField; + CoarseVector Cin(Coarse5D); + CoarseVector Cout(Coarse5D); + CoarseVector CFout(Coarse5D); + + FineField Fin(FineGrid); + FineField Fout(FineGrid); + + + std::vector seeds({1,2,3,4,5}); + GridParallelRNG RNG(Coarse5D); RNG.SeedFixedIntegers(seeds); + + gaussian(RNG,Cin); + PromoteFromSubspace(_Aggregates,Cin,Fin); + ProjectToSubspace(_Aggregates,Cin,Fin); + + std::cout << GridLogMessage<< "************ "<M(Cin,Cout); + this->Project(Cout); + std::cout << GridLogMessage<< " Cout "<Mdag(Cin,Cout); + this->Project(Cout); + std::cout << GridLogMessage<< " Cout "< Directions(void) { return geom.directions;}; + virtual std::vector Displacements(void){ return geom.displacements;}; +}; + +template class SchurSolverWrapper : public LinearFunction { +private: + CheckerBoardedSparseMatrixBase & _Matrix; + SchurRedBlackBase & _Solver; +public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations trick + ///////////////////////////////////////////////////// + SchurSolverWrapper(CheckerBoardedSparseMatrixBase &Matrix, + SchurRedBlackBase &Solver) + : _Matrix(Matrix), _Solver(Solver) {}; + + void operator() (const Field &in, Field &out){ + + _Solver(_Matrix,in,out); // Mdag M out = Mdag in + + } +}; + +template class SolverWrapper : public LinearFunction { +private: + LinearOperatorBase & _Matrix; + OperatorFunction & _Solver; + LinearFunction & _Guess; +public: + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations trick + ///////////////////////////////////////////////////// + SolverWrapper(LinearOperatorBase &Matrix, + OperatorFunction &Solver, + LinearFunction &Guess) + : _Matrix(Matrix), _Solver(Solver), _Guess(Guess) {}; + + void operator() (const Field &in, Field &out){ + + _Guess(in,out); + _Solver(_Matrix,in,out); // Mdag M out = Mdag in + + } +}; + +// Must use a non-hermitian solver +template +class PVdagMLinearOperator : public LinearOperatorBase { + Matrix &_Mat; + Matrix &_PV; +public: + PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){}; + + virtual std::vector Directions(void) { return _Mat.Directions();}; + virtual std::vector Displacements(void){ return _Mat.Displacements();}; + + 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){ + Field tmp(in.Grid()); + _Mat.M(in,tmp); + _PV.Mdag(tmp,out); + } + void AdjOp (const Field &in, Field &out){ + Field tmp(in.Grid()); + _PV.M(tmp,out); + _Mat.Mdag(in,tmp); + } + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + assert(0); + } + void HermOp(const Field &in, Field &out){ + assert(0); + } +}; + +RealD InverseApproximation(RealD x){ + return 1.0/x; +} + +template class ChebyshevSmoother : public LinearFunction +{ +public: + typedef LinearOperatorBase FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; + + Chebyshev Cheby; + + ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) : + _SmootherOperator(SmootherOperator), + _SmootherMatrix(SmootherMatrix), + Cheby(_lo,_hi,_ord,InverseApproximation) + {}; + + void operator() (const Field &in, Field &out) + { + Field tmp(in.Grid()); + MdagMLinearOperator MdagMOp(_SmootherMatrix); + _SmootherOperator.AdjOp(in,tmp); + Cheby(MdagMOp,tmp,out); + } +}; +template +class MGPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + typedef CoarseCayleyFermion CoarseOperator; + // typedef SparseMatrixBase CoarseOperator; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + + MGPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + auto CoarseGrid = _CoarseOperator.Grid(); + CoarseVector Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + std::cout< +class HDCRPreconditioner : public LinearFunction< Lattice > { +public: + + typedef Aggregation Aggregates; + typedef typename Aggregation::CoarseVector CoarseVector; + typedef typename Aggregation::CoarseMatrix CoarseMatrix; + typedef typename Aggregation::FineField FineField; + typedef LinearOperatorBase FineOperator; + typedef LinearFunction FineSmoother; + //typedef CoarseCayleyFermion CoarseOperator; + typedef SparseMatrixBase CoarseOperator; + + Aggregates & _Aggregates; + FineOperator & _FineOperator; + FineSmoother & _PreSmoother; + FineSmoother & _PostSmoother; + CoarseOperator & _CoarseOperator; + CoarseSolver & _CoarseSolve; + + int level; void Level(int lv) {level = lv; }; + + HDCRPreconditioner(Aggregates &Agg, + FineOperator &Fine, + FineSmoother &PreSmoother, + FineSmoother &PostSmoother, + CoarseOperator &CoarseOperator_, + CoarseSolver &CoarseSolve_) + : _Aggregates(Agg), + _FineOperator(Fine), + _PreSmoother(PreSmoother), + _PostSmoother(PostSmoother), + _CoarseOperator(CoarseOperator_), + _CoarseSolve(CoarseSolve_), + level(1) { } + + virtual void operator()(const FineField &in, FineField & out) + { + auto CoarseGrid = _CoarseOperator.Grid(); + CoarseVector Csrc(CoarseGrid); + CoarseVector g5Csrc(CoarseGrid); + CoarseVector Csol(CoarseGrid); + FineField vec1(in.Grid()); + FineField vec2(in.Grid()); + + std::cout< block ({2,2,2,2}); // 4,2,2,2 gets worse + std::vector blockc ({1,1,1,1}); + const int nbasis= 24; + const int nbasisc= 40; // decrease, not improvement + + auto clatt = GridDefaultLatt(); + for(int d=0;d seeds({1,2,3,4}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds); + GridParallelRNG CRNG(Coarse4d);CRNG.SeedFixedIntegers(seeds); + + LatticeGaugeField Umu(UGrid); +#if 0 + SU3::TepidConfiguration(RNG4,Umu); + RealD M5=1.0; +#else + std::string file("./ckpoint_lat.1000"); + FieldMetaData header; + NerscIO::readConfiguration(Umu,header,file); + RealD M5=1.8; +#endif + + std::cout< Subspace; + typedef CoarsenedMatrix CoarseOperator; + typedef CoarseOperator::CoarseVector CoarseVector; + typedef CoarseOperator::siteVector siteVector; + + std::cout< MdagM_Dw(Dw_null); + + std::cout< WilsonCG(1.0e-10,40000); + LatticeFermion w_src(UGrid); w_src=1.0; + LatticeFermion w_res(UGrid); + WilsonCG(MdagM_Dw,w_src,w_res); + exit(0); + */ + std::cout< Level1Op4; + typedef CoarseCayleyFermion Level1Op5; + Level1Op4 c_Dw (*Coarse4d,0); + NonHermitianLinearOperator LinOpDw(Dw); + c_Dw.CoarsenOperator(UGrid,LinOpDw,Aggregates4D); // contains the M5 from Dw(-M5) + // c_Dw.Test(Aggregates4D,UGrid,LinOpDw); + + std::cout< MdagM_cDwf(c_Dwf); + + std::cout<,nbasisc> Level2Op; + typedef Aggregation,nbasisc> CoarseSubspace; + CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0); + + std::cout< L1Hdwf(c_Dwf); + Level2Op cc_Dwf (*CoarseCoarse5d,*CoarseCoarse5dRB,1); // say it is hermitian + cc_Dwf.CoarsenOperator(Coarse5d,L1Hdwf,CoarseAggregates); + // cc_Dwf.Test(CoarseAggregates,Coarse5d,L1Hdwf); + + typedef Level2Op::CoarseVector CoarseCoarseVector; + + std::cout< CoarseCG(tol,MaxIt); + ConjugateGradient FineCG(tol,MaxIt); + + NonHermitianLinearOperator FineM(Ddwf); + MdagMLinearOperator FineMdagM(Ddwf); // M^\dag M + + NonHermitianLinearOperator CoarseM(c_Dwf); + MdagMLinearOperator CoarseMdagM(c_Dwf); + + NonHermitianLinearOperator CoarseCoarseM(cc_Dwf); + MdagMLinearOperator CoarseCoarseMdagM(cc_Dwf); + + + std::cout< PM; PM(MdagM_Dw,w_src); + std::cout< cPM; cPM(CoarseMdagM,c_src); + + cc_src=1.0; + PowerMethod ccPM; ccPM(CoarseCoarseMdagM,cc_src); + + std::cout< IRLHermOpL2(cc_Dwf); + Chebyshev IRLChebyL2(IRL_lo,IRL_hi,IRL_ord); + FunctionHermOp IRLOpChebyL2(IRLChebyL2,IRLHermOpL2); + PlainHermOp IRLOpL2 (IRLHermOpL2); + ImplicitlyRestartedLanczos IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20); + + int cNconv; + cNm=0; + std::vector eval2(cNm); + std::vector evec2(cNm,CoarseCoarse5d); + cc_src=1.0; + // IRLL2.calc(eval2,evec2,cc_src,cNconv); + + std::vector tols ({0.005,0.001}); + std::vector c_los ({0.1,0.05}); + std::vector c_his ({22.0}); + std::vector f_los ({0.5,0.2}); + std::vector f_his ({60.0}); + std::vector ws ({2,3}); + std::vector c_ords ({32,24}); + std::vector f_ords ({20,16}); + + for(auto w : ws ) { + for(auto tol : tols ) { + for(auto f_ord : f_ords ) { + for(auto c_ord : c_ords ) { + for(auto c_lo : c_los ) { + for(auto c_hi : c_his ) { + for(auto f_lo : f_los ) { + for(auto f_hi : f_his ) { + ZeroGuesser CoarseZeroGuesser; + ZeroGuesser CoarseCoarseZeroGuesser; + ConjugateGradient CoarseCoarseCG(tol,10000); + ZeroGuesser CoarseCoarseGuesser; + SchurRedBlackDiagMooeeSolve CoarseCoarseRBCG(CoarseCoarseCG); + SchurSolverWrapper CoarseCoarseSolver(cc_Dwf,CoarseCoarseRBCG); + + std::cout< CoarseCoarseCGNE(cc_Dwf,CoarseCoarseCG,CoarseCoarseZeroGuesser); + { +typedef HDCRPreconditioner,nbasisc,LinearFunction > CoarseMG; + typedef MGPreconditioner > ThreeLevelMG; + + // MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space + // ChebyshevSmoother CoarseSmoother1(0.5,22.0,c_ord,CoarseM,c_Dwf); // 37s, 26 iter + // ChebyshevSmoother CoarseSmoother2(0.5,22.0,c_ord,CoarseM,c_Dwf); + ChebyshevSmoother CoarseSmoother(c_lo,c_hi,c_ord,CoarseM,c_Dwf); // 37s, 26 iter + + // ChebyshevSmoother CoarseSmoother1(0.5,22.0,7,CoarseM,c_Dwf); // 38s, 26 iter + // ChebyshevSmoother CoarseSmoother2(0.5,22.0,7,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.4,22.0,7,CoarseM,c_Dwf); // 41s, 27 iter + // ChebyshevSmoother CoarseSmoother2(0.4,22.0,7,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.6,22.0,6,CoarseM,c_Dwf); // 26 iter + // ChebyshevSmoother CoarseSmoother2(0.6,22.0,6,CoarseM,c_Dwf); + // ChebyshevSmoother CoarseSmoother1(0.5,22.0,5,CoarseM,c_Dwf); // 33 iter, 55s + // ChebyshevSmoother CoarseSmoother2(0.5,22.0,5,CoarseM,c_Dwf); + + + CoarseMG Level2Precon (CoarseAggregates, + CoarseM, + CoarseSmoother, + CoarseSmoother, + cc_Dwf, + CoarseCoarseSolver); + Level2Precon.Level(2); + + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.5, 100, CoarseM,Level2Precon,16,16); // 26 iter, 37s + // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.0, 1, CoarseM,Level2Precon,2,2); // 296 s, 50 iter + // PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.0, 1, CoarseM,Level2Precon,2,2); // 250 s, 37 iter + PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.0, 1, CoarseM,Level2Precon,2,2); + + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(1.0, 100, CoarseM,Level2Precon,16,16); // 35 iter, 45s + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.6, 100, CoarseM,Level2Precon,16,16); // 26,38 (diifferene is measurement noise) + //PrecGeneralisedConjugateResidualNonHermitian L2PGCR(0.2, 100, CoarseM,Level2Precon,16,16); // 26 iter, 47s + L2PGCR.Level(2); + + // Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space + + // ChebyshevSmoother FineSmoother1(0.5,60.0,14,FineM,Ddwf); // 26 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,14,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 25 iter, 38s + // ChebyshevSmoother FineSmoother2(0.5,60.0,16,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 23 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,20,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,10,FineM,Ddwf);24 iter, 44s + // ChebyshevSmoother FineSmoother2(0.5,60.0,24,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // odd convergence tail at 10^-9 ish + // ChebyshevSmoother FineSmoother2(0.1,60.0,24,FineM,Ddwf); // 33 iter, waas O(10-9 by 26) + + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 25 iter, 39s + // ChebyshevSmoother FineSmoother2(0.5,60.0,18,FineM,Ddwf); // + + ChebyshevSmoother FineSmoother(f_lo,f_hi,f_ord,FineM,Ddwf); + + // ChebyshevSmoother FineSmoother1(0.5,60.0,11,FineM,Ddwf); // 33 iter, 49s + // ChebyshevSmoother FineSmoother2(0.5,60.0,11,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.5,60.0,12,FineM,Ddwf); // 26 iter, 37s + // ChebyshevSmoother FineSmoother2(0.5,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.4,60.0,12,FineM,Ddwf); // iter 26 no change in final residual + // ChebyshevSmoother FineSmoother2(0.4,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.3,60.0,12,FineM,Ddwf); // 27 iter 39s. + // ChebyshevSmoother FineSmoother2(0.3,60.0,12,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(0.3,60.0,13,FineM,Ddwf); // 26 iter, but slower + // ChebyshevSmoother FineSmoother2(0.3,60.0,13,FineM,Ddwf); + // ChebyshevSmoother FineSmoother1(1.0,60.0,12,FineM,Ddwf); // 34 iter, slower + // ChebyshevSmoother FineSmoother2(1.0,60.0,12,FineM,Ddwf); + + ThreeLevelMG ThreeLevelPrecon(Aggregates4D, + FineM, + FineSmoother, + FineSmoother, + c_Dwf, + L2PGCR); + ThreeLevelPrecon.Level(1); + + PrecGeneralisedConjugateResidualNonHermitian L1PGCR(1.0e-8,1000,FineM,ThreeLevelPrecon,16,16); + L1PGCR.Level(1); + + f_res=Zero(); + L1PGCR(f_src,f_res); + } + }}}} + }}} + } + std::cout<