From b8a38f292d8f533dedb69cac33cb313f4c31c475 Mon Sep 17 00:00:00 2001 From: paboyle Date: Sat, 28 Nov 2015 17:01:51 -0800 Subject: [PATCH] Domain decomposition SAP precon implemented and working but not as fast as I hoped. --- tests/Test_dwf_hdcr.cc | 317 ++++++++++++++++++++++++++++++----------- 1 file changed, 231 insertions(+), 86 deletions(-) diff --git a/tests/Test_dwf_hdcr.cc b/tests/Test_dwf_hdcr.cc index 94b8481a..eb247331 100644 --- a/tests/Test_dwf_hdcr.cc +++ b/tests/Test_dwf_hdcr.cc @@ -6,6 +6,22 @@ using namespace std; using namespace Grid; using namespace Grid::QCD; +class myclass { +public: + + GRID_DECL_CLASS_MEMBERS(myclass, + int, domaindecompose, + int, domainsize, + int, order, + double, lo, + double, hi, + int, steps); + + myclass(){}; + +}; +myclass params; + RealD InverseApproximation(RealD x){ return 1.0/x; } @@ -26,15 +42,21 @@ public: Aggregates & _Aggregates; CoarseOperator & _CoarseOperator; - Matrix & _Matrix; + Matrix & _FineMatrix; FineOperator & _FineOperator; + Matrix & _SmootherMatrix; + FineOperator & _SmootherOperator; // Constructor - MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, FineOperator &Fine,Matrix &FineMatrix) + MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, + FineOperator &Fine,Matrix &FineMatrix, + FineOperator &Smooth,Matrix &SmootherMatrix) : _Aggregates(Agg), _CoarseOperator(Coarse), _FineOperator(Fine), - _Matrix(FineMatrix) + _FineMatrix(FineMatrix), + _SmootherOperator(Smooth), + _SmootherMatrix(SmootherMatrix) { } @@ -43,7 +65,7 @@ public: FineField p1(in._grid); FineField p2(in._grid); - MdagMLinearOperator fMdagMOp(_Matrix); + MdagMLinearOperator fMdagMOp(_FineMatrix); p1=in; RealD absp2; @@ -58,74 +80,20 @@ public: } } -#if 0 void operator()(const FineField &in, FineField & out) { - - FineField Min(in._grid); - FineField tmp(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< fCG(1.0e-3,1000); - ConjugateGradient CG(1.0e-8,100000); - - //////////////////////////////////////////////////////////////////////// - // ADEF2: [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] - //////////////////////////////////////////////////////////////////////// - - // Smoothing step, followed by coarse grid correction - MdagMLinearOperator MdagMOp(_Matrix); - - Min=in; - std::cout< 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;; - */ - + if ( params.domaindecompose ) { + operatorSAP(in,out); + } else { + operatorCheby(in,out); + } } -#endif //////////////////////////////////////////////////////////////////////// // ADEF2: [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] // ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in //////////////////////////////////////////////////////////////////////// -#if 0 - void operator()(const FineField &in, FineField & out) { +#if 1 + void operatorADEF2(const FineField &in, FineField & out) { CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid()); @@ -136,7 +104,7 @@ public: HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); - MdagMLinearOperator fMdagMOp(_Matrix); + MdagMLinearOperator fMdagMOp(_FineMatrix); FineField tmp(in._grid); FineField res(in._grid); @@ -189,8 +157,8 @@ public: } #endif // ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in -#if 0 - void operator()(const FineField &in, FineField & out) { +#if 1 + void operatorADEF1(const FineField &in, FineField & out) { CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid()); @@ -201,7 +169,7 @@ public: HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); - ShiftedMdagMLinearOperator fMdagMOp(_Matrix,0.1); + ShiftedMdagMLinearOperator fMdagMOp(_FineMatrix,0.1); FineField tmp(in._grid); FineField res(in._grid); @@ -234,14 +202,79 @@ public: } #endif + void SAP (const FineField & src,FineField & psi){ + + Lattice > coor(src._grid); + Lattice > subset(src._grid); + + FineField r(src._grid); + FineField zz(src._grid); zz=zero; + FineField vec1(src._grid); + FineField vec2(src._grid); + + const Integer block=params.domainsize; + + subset=zero; + for(int mu=0;mu fMdagMOp(_SmootherMatrix,0.0); + Chebyshev Cheby (params.lo,params.hi,params.order,InverseApproximation); + + RealD resid; + for(int i=0;i 20*71 = 1400 matmuls. +// 2*71 = 140 comms. + + // Even domain solve + r= where(subset==(Integer)0,r,zz); + _SmootherOperator.AdjOp(r,vec1); + Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M + psi = psi + vec2; + + // Odd domain residual + _FineOperator.Op(psi,vec1);// this is the G5 herm bit + r= src - vec1 ; + r= where(subset==(Integer)1,r,zz); + + resid = norm2(r) /norm2(src); + std::cout << "SAP "< fMdagMOp(_Matrix); - ShiftedMdagMLinearOperator fMdagMOp(_Matrix,0.5); + // MdagMLinearOperator fMdagMOp(_FineMatrix); + ShiftedMdagMLinearOperator fMdagMOp(_SmootherMatrix,0.0); RealD Ni,r; @@ -250,7 +283,7 @@ public: for(int ilo=0;ilo<3;ilo++){ for(int ord=5;ord<50;ord*=2){ - _FineOperator.AdjOp(in,vec1); + _SmootherOperator.AdjOp(in,vec1); Chebyshev Cheby (lo[ilo],70.0,ord,InverseApproximation); Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M @@ -264,7 +297,7 @@ public: } } - void operator()(const FineField &in, FineField & out) { + void operatorCheby(const FineField &in, FineField & out) { CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid()); @@ -275,18 +308,18 @@ public: HermitianLinearOperator HermOp(_CoarseOperator); MdagMLinearOperator MdagMOp(_CoarseOperator); - // MdagMLinearOperator fMdagMOp(_Matrix); - ShiftedMdagMLinearOperator fMdagMOp(_Matrix,0.0); + // MdagMLinearOperator fMdagMOp(_FineMatrix); + ShiftedMdagMLinearOperator fMdagMOp(_SmootherMatrix,0.0); FineField vec1(in._grid); FineField vec2(in._grid); // Chebyshev Cheby (0.5,70.0,30,InverseApproximation); // Chebyshev ChebyAccu(0.5,70.0,30,InverseApproximation); - Chebyshev Cheby (2.0,70.0,10,InverseApproximation); - Chebyshev ChebyAccu(2.0,70.0,10,InverseApproximation); - Cheby.JacksonSmooth(); - ChebyAccu.JacksonSmooth(); + Chebyshev Cheby (2.0,70.0,15,InverseApproximation); + Chebyshev ChebyAccu(2.0,70.0,15,InverseApproximation); + // Cheby.JacksonSmooth(); + // ChebyAccu.JacksonSmooth(); _Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.PromoteFromSubspace(Csrc,out); @@ -305,7 +338,7 @@ public: RealD Ni = norm2(in); - _FineOperator.AdjOp(in,vec1);// this is the G5 herm bit + _SmootherOperator.AdjOp(in,vec1);// this is the G5 herm bit ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M std::cout< CG(1.0e-3,100000); + + HermitianLinearOperator HermOp(_CoarseOperator); + MdagMLinearOperator MdagMOp(_CoarseOperator); + + FineField vec1(in._grid); + FineField vec2(in._grid); + + _Aggregates.ProjectToSubspace (Csrc,in); + _Aggregates.PromoteFromSubspace(Csrc,out); + std::cout< > 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); @@ -402,6 +517,7 @@ int main (int argc, char ** argv) std::cout< HermIndefOp(Ddwf); + Gamma5R5HermitianLinearOperator HermIndefOpDD(DdwfDD); CoarsenedMatrix LDOp(*Coarse5d); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); @@ -467,7 +584,13 @@ int main (int argc, char ** argv) std::cout< Precon(Aggregates, LDOp,HermIndefOp,Ddwf); + MultiGridPreconditioner Precon (Aggregates, LDOp, + HermIndefOp,Ddwf, + HermIndefOp,Ddwf); + + MultiGridPreconditioner PreconDD(Aggregates, LDOp, + HermIndefOp,Ddwf, + HermIndefOpDD,DdwfDD); TrivialPrecon simple; std::cout< simple; // ConjugateGradient fCG(1.0e-8,100000); // fCG(HermDefOp,src,result); @@ -496,12 +630,22 @@ int main (int argc, char ** argv) std::cout< PGCRDD(1.0e-8,100000,PreconDD,8,128); + result=zero; + std::cout< PGCR(1.0e-8,100000,Precon,8,128); - std::cout< PGCR(1.0e-8,100000,Precon,8,128); + // std::cout<