diff --git a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h index af043ff4c..aef53cc72 100644 --- a/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h +++ b/Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h @@ -28,6 +28,7 @@ Author: Peter Boyle /* END LEGAL */ #ifndef GRID_PREC_GCR_NON_HERM_H #define GRID_PREC_GCR_NON_HERM_H +#include /////////////////////////////////////////////////////////////////////////////////////////////////////// //VPGCR Abe and Zhang, 2005. @@ -79,6 +80,11 @@ public: // same applies and no reductions. Off by default; boss rank prints. int LogCoeffs = 0; void LogCoefficients(int l) { LogCoeffs = l; }; + // Optional recorder of the per-step coefficients (means over calls), for + // replay by GCRReplaySmoother (Smoothers.h). Records only; no effect on + // the iteration. + GCRCoefficients *Recorder = nullptr; + void SetCoefficientRecorder(GCRCoefficients *r) { Recorder = r; if(r) r->mmax = mmax; }; PrecGeneralisedConjugateResidualNonHermitian(RealD tol,Integer maxit,LinearOperatorBase &_Linop,LinearFunction &Prec,int _mmax,int _nstep) : Tolerance(tol), @@ -223,6 +229,7 @@ public: LinalgTimer.Start(); rq= innerProduct(q[peri_k],r); // what if rAr not real? a = rq/qq[peri_k]; + if ( Recorder ) Recorder->RecordA(k,a); axpy(psi,a,p[peri_k],psi); @@ -274,6 +281,7 @@ public: bcoef[back] = -bcoef[back]/qq[peri_back]; if ( LogCoeffs ) bs<<" b["<RecordB(k,bcoef); if ( northog ) { axpyMulti(p[peri_kp],bcoef,pwin); qq[peri_kp]=axpyMultiNorm(q[peri_kp],bcoef,qwin); diff --git a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h index 20556c346..7ae154623 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h +++ b/Grid/algorithms/multigrid/BlockCyclicSchurInverse.h @@ -140,6 +140,7 @@ public: /////////////////////////////////////////////////////////////////////////// void Leaf(BlockCyclicMatrix &A, int64_t b) { + GRID_TRACE("SchurLeaf"); BlockCyclicLayout &L = A.layout; nLeaf++; if ( (int)(b % L.Pr) != L.prow ) return; @@ -211,6 +212,7 @@ public: int64_t span = b1-b0; GRID_ASSERT( span >= 1 ); if ( span == 1 ) { Leaf(A, b0); return; } + GRID_TRACE("SchurNode"); nNode++; int64_t bm = b0 + span/2; @@ -244,8 +246,10 @@ public: tGemm += usecond(); // 9. Off-diagonal signs + { GRID_TRACE("SchurCopy"); WindowCopyScale(mone, Ut, A, c0,m, m,c1); WindowCopyScale(mone, Tt, A, m,c1, c0,m); + } } /////////////////////////////////////////////////////////////////////////// @@ -285,6 +289,22 @@ public: << " leaf " << tLeaf/1.0e6 << " copy " << tCopy/1.0e6 << ")" << std::endl; + // SUMMA breakdown: boss-rank seconds, plus the ring/gemm time spread over + // ranks (min/max) -- skew shows as max >> min. + RealD ring = (SUMMA.tRingA+SUMMA.tRingB)/1.0e6; + RealD rmin = -ring, rmax = ring; grid->GlobalMax(rmin); grid->GlobalMax(rmax); rmin = -rmin; // no GlobalMin: max of negation + RealD gmin = -SUMMA.tGemm/1.0e6, gmax = SUMMA.tGemm/1.0e6; grid->GlobalMax(gmin); grid->GlobalMax(gmax); gmin = -gmin; + double gb = SUMMA.bytesRing/1.0e9; + std::cout << GridLogMessage << "BlockCyclicSumma:" + << " multiplies " << SUMMA.nMultiply << " gemms " << SUMMA.nGemm << " ring msgs " << SUMMA.nRingMsg + << " | boss secs: alloc " << SUMMA.tAlloc/1.0e6 + << " pack " << SUMMA.tPack/1.0e6 + << " ringA " << SUMMA.tRingA/1.0e6 << " ringB " << SUMMA.tRingB/1.0e6 + << " gemm " << SUMMA.tGemm/1.0e6 + << " | ring min/max over ranks " << rmin << "/" << rmax + << " gemm min/max " << gmin << "/" << gmax + << " | ring bytes/rank " << gb << " GB -> " << (ring>0 ? gb/ring : 0.0) << " GB/s/rank (boss)" + << std::endl; } }; diff --git a/Grid/algorithms/multigrid/BlockCyclicSumma.h b/Grid/algorithms/multigrid/BlockCyclicSumma.h index 2359afcd6..625667e58 100644 --- a/Grid/algorithms/multigrid/BlockCyclicSumma.h +++ b/Grid/algorithms/multigrid/BlockCyclicSumma.h @@ -127,6 +127,15 @@ class BlockCyclicSumma public: GridBLAS BLAS; + // Per-rank telemetry (boss-rank seconds when printed; no reductions here). + // Every Multiply is: buffer alloc, pack panels, ring A along the process + // row, ring B along the process column, then the local GEMMs. The rings + // are synchronous SendToRecvFrom, so tRing is time the GPU is idle unless + // a future version overlaps them with the GEMMs. + double tAlloc=0, tPack=0, tRingA=0, tRingB=0, tGemm=0; + uint64_t bytesRing=0, nRingMsg=0, nMultiply=0, nGemm=0; + void ResetTelemetry(void){ tAlloc=tPack=tRingA=tRingB=tGemm=0; bytesRing=nRingMsg=nMultiply=nGemm=0; } + static int Overlap(int64_t a0,int64_t a1,int64_t b0,int64_t b1) { return (a0 < b1) && (b0 < a1); } @@ -188,8 +197,11 @@ public: const uint64_t slotA = (uint64_t)mloc_i*nb; const uint64_t slotB1 = (uint64_t)nb*nloc_j; // one panel const uint64_t slotB = (uint64_t)S*slotB1; + nMultiply++; + tAlloc -= usecond(); deviceVector Abuf( slotA*Pc ? slotA*Pc : 1 ); deviceVector Bbuf( slotB*Pr ? slotB*Pr : 1 ); + tAlloc += usecond(); deviceVector ap(1), bp(1), cp(1); std::vector ptr(1); @@ -201,6 +213,8 @@ public: ///////////////////////////////////////////////////////////////////// // Pack MY panels of this round into my origin slots. ///////////////////////////////////////////////////////////////////// + tPack -= usecond(); + { GRID_TRACE("SummaPack"); for(int64_t s=r0; s 1 && slotA ){ + GRID_TRACE("SummaRingA"); + tRingA -= usecond(); int dest = prow*Pc + (pcol+1)%Pc; int src = prow*Pc + (pcol-1+Pc)%Pc; for(int t=1;tSendToRecvFrom((void *)(&Abuf[0]+slotA*cs), dest, (void *)(&Abuf[0]+slotA*cr), src, slotA*sizeof(ComplexD)); + bytesRing += slotA*sizeof(ComplexD); nRingMsg++; } + tRingA += usecond(); } ///////////////////////////////////////////////////////////////////// // Ring allgather along my process COLUMN: Pr-1 symmetric steps. ///////////////////////////////////////////////////////////////////// if ( Pr > 1 && slotB ){ + GRID_TRACE("SummaRingB"); + tRingB -= usecond(); int dest = ((prow+1)%Pr)*Pc + pcol; int src = ((prow-1+Pr)%Pr)*Pc + pcol; for(int t=1;tSendToRecvFrom((void *)(&Bbuf[0]+slotB*rs), dest, (void *)(&Bbuf[0]+slotB*rr), src, slotB*sizeof(ComplexD)); + bytesRing += slotB*sizeof(ComplexD); nRingMsg++; } + tRingB += usecond(); } ///////////////////////////////////////////////////////////////////// // Local update, ascending s: fixed order, bitwise-reproducible. ///////////////////////////////////////////////////////////////////// + tGemm -= usecond(); + { GRID_TRACE("SummaGEMM"); for(int64_t s=r0; s #include #include #include +#include diff --git a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc index 18fbe38a6..c38f6cd7d 100644 --- a/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc +++ b/examples/Example_pvdagm_v2_3level_DenseCoarseMatrix.cc @@ -93,6 +93,18 @@ int FineSmootherOrder = 6; int FineSmootherMmax = 6; RealD CoarseSmootherShift = 0.1; int PowerIterations = 0; // >0: power-iterate the smoother operators before the solves (spectral edge) +// Smoother implementation per level (Smoothers.h): +// gcr : the adaptive PGCR (default, as always) +// replay : run the PGCR with a coefficient recorder for the first +// PolyRecordIters outer steps, then switch to GCRReplaySmoother +// (same polynomial, no inner products) -- 1402.2585 p.13 revisited +// cheb : ChebyshevNonHermitianSmoother, 1/x on [ChebLo,ChebHi], order +// = the GCR step count of that level +std::string FineSmootherMode = "gcr"; +std::string CoarseSmootherMode = "gcr"; +int PolyRecordIters = 4; +RealD FineChebLo = 3.0, FineChebHi = 137.0; // from the harvested polynomial and PowerIterations edge +RealD CoarseChebLo = 8.0, CoarseChebHi = 45.0; int CoarseSmootherNstep = 2; int CoarseSmootherMmax = 2; RealD CoarseSolverTol = 0.05; @@ -137,6 +149,13 @@ void ParseEnvironment(void) if(getenv("FineSmootherMmax")) FineSmootherMmax = atoi(getenv("FineSmootherMmax")); if(getenv("CoarseSmootherShift"))CoarseSmootherShift= atof(getenv("CoarseSmootherShift")); if(getenv("PowerIterations")) PowerIterations = atoi(getenv("PowerIterations")); + if(getenv("FineSmootherMode")) FineSmootherMode = getenv("FineSmootherMode"); + if(getenv("CoarseSmootherMode")) CoarseSmootherMode = getenv("CoarseSmootherMode"); + if(getenv("PolyRecordIters")) PolyRecordIters = atoi(getenv("PolyRecordIters")); + if(getenv("FineChebLo")) FineChebLo = atof(getenv("FineChebLo")); + if(getenv("FineChebHi")) FineChebHi = atof(getenv("FineChebHi")); + if(getenv("CoarseChebLo")) CoarseChebLo = atof(getenv("CoarseChebLo")); + if(getenv("CoarseChebHi")) CoarseChebHi = atof(getenv("CoarseChebHi")); if(getenv("CoarseSmootherNstep"))CoarseSmootherNstep= atoi(getenv("CoarseSmootherNstep")); if(getenv("CoarseSmootherMmax")) CoarseSmootherMmax = atoi(getenv("CoarseSmootherMmax")); if(getenv("CoarseSolverTol")) CoarseSolverTol = atof(getenv("CoarseSolverTol")); @@ -346,6 +365,7 @@ public: std::string name = "Level 1"; LinearOperatorBase &Linop; MrhsLinearFunction &Preconditioner; + std::function OnStep; // called with the outer step count after every step (smoother switching) void Level(int lv){ name = "Level " + std::to_string(lv); level=lv; } void Name(std::string n){ name = n; } void SetZeroGuess(int z){ ZeroGuess=z; } @@ -387,6 +407,7 @@ public: Preconditioner(r,p[0]); vOp(p[0],q[0]); qq[0]=vnorm2(q[0]); cp=vnorm2(r); for(int k=0;k CoarseSmootherSlot(CoarseSmootherGCR,"Csmoother GCR"); MrhsCoarseThreeLevelPrec - L2to3Precon(LinOpC, CoarseSmootherGCR, MrhsProjectorL2, ccSolve, + L2to3Precon(LinOpC, CoarseSmootherSlot, MrhsProjectorL2, ccSolve, Coarse5d, CoarseCoarse5d, CCMrhs, nr); PrecGeneralisedConjugateResidualNonHermitian @@ -951,13 +973,51 @@ int main (int argc, char ** argv) SmootherGCR.LogCoefficients(SmootherCoeffLog); CoarseSmootherGCR.LogCoefficients(SmootherCoeffLog); - MrhsTwoLevelMG - ThreeLevelPrecon(PVdagM, SmootherGCR, MrhsProjector, L2PGCR, Coarse5d, CMrhs); + SwitchableSmoother FineSmootherSlot(SmootherGCR,"Fsmoother GCR"); + MrhsTwoLevelMG > + ThreeLevelPrecon(PVdagM, FineSmootherSlot, MrhsProjector, L2PGCR, Coarse5d, CMrhs); MrhsPGCRNonHermitian L1PGCR(OuterTol,1000,PVdagM,ThreeLevelPrecon,OuterMmax,OuterNstep); L1PGCR.Level(1); L1PGCR.Name("Fouter"); L1PGCR.SetZeroGuess(1); + ////////////////////////////////////////////////////////////////////// + // Smoother modes. Objects live for the duration of this RunSolve. + ////////////////////////////////////////////////////////////////////// + std::unique_ptr > FineCheb; + std::unique_ptr > CoarseCheb; + std::unique_ptr > FineReplay; + std::unique_ptr > CoarseReplay; + GCRCoefficients recF, recC; + if ( FineSmootherMode == "cheb" ) { + FineCheb.reset(new ChebyshevNonHermitianSmoother(FineChebLo,FineChebHi,FineSmootherOrder,ShiftedPVdagM)); + FineSmootherSlot.Set(*FineCheb,"Fsmoother Chebyshev"); + } + if ( CoarseSmootherMode == "cheb" ) { + CoarseCheb.reset(new ChebyshevNonHermitianSmoother(CoarseChebLo,CoarseChebHi,CoarseSmootherNstep,ShiftedC)); + CoarseSmootherSlot.Set(*CoarseCheb,"Csmoother Chebyshev"); + } + if ( FineSmootherMode == "replay" ) SmootherGCR.SetCoefficientRecorder(&recF); + if ( CoarseSmootherMode == "replay" ) CoarseSmootherGCR.SetCoefficientRecorder(&recC); + L1PGCR.OnStep = [&](int step){ + if ( step != PolyRecordIters ) return; + if ( FineSmootherMode == "replay" ) { + SmootherGCR.SetCoefficientRecorder(nullptr); + recF.Report("Fsmoother"); + FineReplay.reset(new GCRReplaySmoother(ShiftedPVdagM,recF)); + FineSmootherSlot.Set(*FineReplay,"Fsmoother replay"); + } + if ( CoarseSmootherMode == "replay" ) { + CoarseSmootherGCR.SetCoefficientRecorder(nullptr); + recC.Report("Csmoother"); + CoarseReplay.reset(new GCRReplaySmoother(ShiftedC,recC)); + CoarseSmootherSlot.Set(*CoarseReplay,"Csmoother replay"); + } + }; + std::cout << GridLogMessage << "Smoother modes: fine " << FineSmootherMode << " coarse " << CoarseSmootherMode + << (FineSmootherMode=="replay"||CoarseSmootherMode=="replay" ? " (record for "+std::to_string(PolyRecordIters)+" outer steps)" : "") + << std::endl; + std::vector src(nr,FGrid), sol(nr,FGrid); for(int r=0;r @@ -132,19 +129,6 @@ public: void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ GRID_ASSERT(0); } }; /* -template class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - Chebyshev Cheby; - ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - Cheby(_lo,_hi,_ord,InverseApproximation) - { - std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - Chebyshev Cheby; - ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - Cheby(_lo,_hi,_ord,InverseApproximation) - { - std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class ChebyshevInverter : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _Operator; - Chebyshev Cheby; - ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) : - _Operator(Operator), - Cheby(_lo,_hi,_ord,InverseApproximation) - { - std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _SmootherOperator; - Chebyshev Cheby; - ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) : - _SmootherOperator(SmootherOperator), - Cheby(_lo,_hi,_ord,InverseApproximation) - { - std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class ChebyshevInverter : public LinearFunction -{ -public: - using LinearFunction::operator(); - typedef LinearOperatorBase FineOperator; - FineOperator & _Operator; - Chebyshev Cheby; - ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) : - _Operator(Operator), - Cheby(_lo,_hi,_ord,InverseApproximation) - { - std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"< class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 > { diff --git a/tests/solver/Test_dwf_hdcr_16_rb.cc b/tests/solver/Test_dwf_hdcr_16_rb.cc index ee613aadc..2ae0d2e48 100644 --- a/tests/solver/Test_dwf_hdcr_16_rb.cc +++ b/tests/solver/Test_dwf_hdcr_16_rb.cc @@ -48,9 +48,6 @@ using namespace Grid; * 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: @@ -72,60 +69,6 @@ public: } }; -template class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 > { diff --git a/tests/solver/Test_dwf_hdcr_24_regression.cc b/tests/solver/Test_dwf_hdcr_24_regression.cc index 023a9e8b0..c28282983 100644 --- a/tests/solver/Test_dwf_hdcr_24_regression.cc +++ b/tests/solver/Test_dwf_hdcr_24_regression.cc @@ -48,64 +48,7 @@ using namespace Grid; * 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: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 > { diff --git a/tests/solver/Test_dwf_hdcr_2level.cc b/tests/solver/Test_dwf_hdcr_2level.cc index 0975fa402..28168463d 100644 --- a/tests/solver/Test_dwf_hdcr_2level.cc +++ b/tests/solver/Test_dwf_hdcr_2level.cc @@ -49,64 +49,7 @@ using namespace Grid; * 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: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 RedBlackSmoother : public LinearFunction { public: diff --git a/tests/solver/Test_dwf_hdcr_48_rb.cc b/tests/solver/Test_dwf_hdcr_48_rb.cc index 3fef75ca8..5de718d1b 100644 --- a/tests/solver/Test_dwf_hdcr_48_rb.cc +++ b/tests/solver/Test_dwf_hdcr_48_rb.cc @@ -48,9 +48,6 @@ using namespace Grid; * 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: @@ -72,60 +69,6 @@ public: } }; -template class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 > { diff --git a/tests/solver/Test_dwf_hdcr_48_regression.cc b/tests/solver/Test_dwf_hdcr_48_regression.cc index 616ab7ca6..54ed879c7 100644 --- a/tests/solver/Test_dwf_hdcr_48_regression.cc +++ b/tests/solver/Test_dwf_hdcr_48_regression.cc @@ -48,64 +48,7 @@ using namespace Grid; * 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: - using LinearFunction::operator(); - 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: - using LinearFunction::operator(); - 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 > { diff --git a/tests/solver/Test_dwf_multigrid.cc b/tests/solver/Test_dwf_multigrid.cc index 0f933971c..bed661ace 100644 --- a/tests/solver/Test_dwf_multigrid.cc +++ b/tests/solver/Test_dwf_multigrid.cc @@ -112,64 +112,8 @@ public: }; -RealD InverseApproximation(RealD x){ - return 1.0/x; -} -template class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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); - } -}; #define GridLogLevel std::cout << GridLogMessage < 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); - } -}; #define GridLogLevel std::cout << GridLogMessage < class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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: diff --git a/tests/solver/Test_hw_multigrid_mixed_48_rb.cc b/tests/solver/Test_hw_multigrid_mixed_48_rb.cc index d8a06f807..c889ae89f 100644 --- a/tests/solver/Test_hw_multigrid_mixed_48_rb.cc +++ b/tests/solver/Test_hw_multigrid_mixed_48_rb.cc @@ -725,34 +725,7 @@ public: } }; -RealD InverseApproximation(RealD x){ - return 1.0/x; -} -template class ChebyshevSmoother : public LinearFunction -{ -public: - using LinearFunction::operator(); - 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: