diff --git a/lib/algorithms/iterative/ConjugateGradientMixedPrec.h b/lib/algorithms/iterative/ConjugateGradientMixedPrec.h index c7332455..0c24a837 100644 --- a/lib/algorithms/iterative/ConjugateGradientMixedPrec.h +++ b/lib/algorithms/iterative/ConjugateGradientMixedPrec.h @@ -1,4 +1,4 @@ - /************************************************************************************* +/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid @@ -23,132 +23,134 @@ Author: Christopher Kelly 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 */ +*************************************************************************************/ +/* END LEGAL */ #ifndef GRID_CONJUGATE_GRADIENT_MIXED_PREC_H #define GRID_CONJUGATE_GRADIENT_MIXED_PREC_H -namespace Grid { +NAMESPACE_BEGIN(Grid); - //Mixed precision restarted defect correction CG - template::value == 2, int>::type = 0,typename std::enable_if< getPrecision::value == 1, int>::type = 0> - class MixedPrecisionConjugateGradient : public LinearFunction { - public: - RealD Tolerance; - RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed - Integer MaxInnerIterations; - Integer MaxOuterIterations; - GridBase* SinglePrecGrid; //Grid for single-precision fields - RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance - LinearOperatorBase &Linop_f; - LinearOperatorBase &Linop_d; +//Mixed precision restarted defect correction CG +template::value == 2, int>::type = 0, + typename std::enable_if< getPrecision::value == 1, int>::type = 0> +class MixedPrecisionConjugateGradient : public LinearFunction { +public: + RealD Tolerance; + RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed + Integer MaxInnerIterations; + Integer MaxOuterIterations; + GridBase* SinglePrecGrid; //Grid for single-precision fields + RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance + LinearOperatorBase &Linop_f; + LinearOperatorBase &Linop_d; - Integer TotalInnerIterations; //Number of inner CG iterations - Integer TotalOuterIterations; //Number of restarts - Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step + Integer TotalInnerIterations; //Number of inner CG iterations + Integer TotalOuterIterations; //Number of restarts + Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step - //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess - LinearFunction *guesser; + //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess + LinearFunction *guesser; - MixedPrecisionConjugateGradient(RealD tol, Integer maxinnerit, Integer maxouterit, GridBase* _sp_grid, LinearOperatorBase &_Linop_f, LinearOperatorBase &_Linop_d) : - Linop_f(_Linop_f), Linop_d(_Linop_d), - Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid), - OuterLoopNormMult(100.), guesser(NULL){ }; + MixedPrecisionConjugateGradient(RealD tol, Integer maxinnerit, Integer maxouterit, GridBase* _sp_grid, LinearOperatorBase &_Linop_f, LinearOperatorBase &_Linop_d) : + Linop_f(_Linop_f), Linop_d(_Linop_d), + Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid), + OuterLoopNormMult(100.), guesser(NULL){ }; - void useGuesser(LinearFunction &g){ - guesser = &g; - } + void useGuesser(LinearFunction &g){ + guesser = &g; + } - void operator() (const FieldD &src_d_in, FieldD &sol_d){ - TotalInnerIterations = 0; + void operator() (const FieldD &src_d_in, FieldD &sol_d){ + TotalInnerIterations = 0; - GridStopWatch TotalTimer; - TotalTimer.Start(); + GridStopWatch TotalTimer; + TotalTimer.Start(); - int cb = src_d_in.checkerboard; - sol_d.checkerboard = cb; + int cb = src_d_in.checkerboard; + sol_d.checkerboard = cb; - RealD src_norm = norm2(src_d_in); - RealD stop = src_norm * Tolerance*Tolerance; + RealD src_norm = norm2(src_d_in); + RealD stop = src_norm * Tolerance*Tolerance; - GridBase* DoublePrecGrid = src_d_in._grid; - FieldD tmp_d(DoublePrecGrid); - tmp_d.checkerboard = cb; + GridBase* DoublePrecGrid = src_d_in._grid; + FieldD tmp_d(DoublePrecGrid); + tmp_d.checkerboard = cb; - FieldD tmp2_d(DoublePrecGrid); - tmp2_d.checkerboard = cb; + FieldD tmp2_d(DoublePrecGrid); + tmp2_d.checkerboard = cb; - FieldD src_d(DoublePrecGrid); - src_d = src_d_in; //source for next inner iteration, computed from residual during operation + FieldD src_d(DoublePrecGrid); + src_d = src_d_in; //source for next inner iteration, computed from residual during operation - RealD inner_tol = InnerTolerance; + RealD inner_tol = InnerTolerance; - FieldF src_f(SinglePrecGrid); - src_f.checkerboard = cb; + FieldF src_f(SinglePrecGrid); + src_f.checkerboard = cb; - FieldF sol_f(SinglePrecGrid); - sol_f.checkerboard = cb; + FieldF sol_f(SinglePrecGrid); + sol_f.checkerboard = cb; - ConjugateGradient CG_f(inner_tol, MaxInnerIterations); - CG_f.ErrorOnNoConverge = false; + ConjugateGradient CG_f(inner_tol, MaxInnerIterations); + CG_f.ErrorOnNoConverge = false; - GridStopWatch InnerCGtimer; + GridStopWatch InnerCGtimer; - GridStopWatch PrecChangeTimer; + GridStopWatch PrecChangeTimer; - Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count + Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count - for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){ - //Compute double precision rsd and also new RHS vector. - Linop_d.HermOp(sol_d, tmp_d); - RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector + for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){ + //Compute double precision rsd and also new RHS vector. + Linop_d.HermOp(sol_d, tmp_d); + RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector - std::cout< CG_d(Tolerance, MaxInnerIterations); - CG_d(Linop_d, src_d_in, sol_d); - TotalFinalStepIterations = CG_d.IterationsToComplete; + while(norm * inner_tol * inner_tol < stop) inner_tol *= 2; // inner_tol = sqrt(stop/norm) ?? - TotalTimer.Stop(); - std::cout< CG_d(Tolerance, MaxInnerIterations); + CG_d(Linop_d, src_d_in, sol_d); + TotalFinalStepIterations = CG_d.IterationsToComplete; -} + TotalTimer.Stop(); + std::cout<