From 0d5af667d87459a24a969979a922ba16d50719d7 Mon Sep 17 00:00:00 2001 From: Chulwoo Jung Date: Tue, 21 Mar 2017 12:40:43 -0400 Subject: [PATCH] Works with CPS evolution --- lib/Algorithms.h | 1 + lib/algorithms/iterative/ConjugateGradient.h | 5 +- .../iterative/ConjugateGradientMixedPrec.h | 32 +- .../iterative/ConjugateGradientMultiShift.h | 11 +- .../ConjugateGradientMultiShiftMixedPrec.h | 404 ++++++++++++++++++ .../iterative/ConjugateGradientShifted.h | 168 ++++++++ .../iterative/ImplicitlyRestartedLanczos.h | 1 + lib/algorithms/iterative/Matrix.h | 2 +- 8 files changed, 597 insertions(+), 27 deletions(-) create mode 100644 lib/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h create mode 100644 lib/algorithms/iterative/ConjugateGradientShifted.h diff --git a/lib/Algorithms.h b/lib/Algorithms.h index 67eb11c3..d586df9d 100644 --- a/lib/Algorithms.h +++ b/lib/Algorithms.h @@ -39,6 +39,7 @@ Author: Peter Boyle #include #include +#include #include #include #include diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h index 0f24ae94..cf3872c8 100644 --- a/lib/algorithms/iterative/ConjugateGradient.h +++ b/lib/algorithms/iterative/ConjugateGradient.h @@ -45,8 +45,6 @@ class ConjugateGradient : public OperatorFunction { // Defaults true. RealD Tolerance; Integer MaxIterations; - Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion - ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true) : Tolerance(tol), MaxIterations(maxit), @@ -157,14 +155,13 @@ class ConjugateGradient : public OperatorFunction { std::cout << std::endl; if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0); - IterationsToComplete = k; + return; } } std::cout << GridLogMessage << "ConjugateGradient did NOT converge" << std::endl; if (ErrorOnNoConverge) assert(0); - IterationsToComplete = k; } }; } diff --git a/lib/algorithms/iterative/ConjugateGradientMixedPrec.h b/lib/algorithms/iterative/ConjugateGradientMixedPrec.h index c7332455..641d58bb 100644 --- a/lib/algorithms/iterative/ConjugateGradientMixedPrec.h +++ b/lib/algorithms/iterative/ConjugateGradientMixedPrec.h @@ -35,7 +35,6 @@ namespace Grid { 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 @@ -43,16 +42,12 @@ namespace Grid { 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 - //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), + Tolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid), OuterLoopNormMult(100.), guesser(NULL){ }; void useGuesser(LinearFunction &g){ @@ -60,8 +55,9 @@ namespace Grid { } void operator() (const FieldD &src_d_in, FieldD &sol_d){ - TotalInnerIterations = 0; - + (*this)(src_d_in,sol_d,NULL); + } + void operator() (const FieldD &src_d_in, FieldD &sol_d, RealD *shift){ GridStopWatch TotalTimer; TotalTimer.Start(); @@ -81,7 +77,7 @@ namespace Grid { 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 = Tolerance; FieldF src_f(SinglePrecGrid); src_f.checkerboard = cb; @@ -89,18 +85,17 @@ namespace Grid { FieldF sol_f(SinglePrecGrid); sol_f.checkerboard = cb; - ConjugateGradient CG_f(inner_tol, MaxInnerIterations); + ConjugateGradientShifted CG_f(inner_tol, MaxInnerIterations); CG_f.ErrorOnNoConverge = false; GridStopWatch InnerCGtimer; GridStopWatch PrecChangeTimer; - Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count - - for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){ + for(Integer outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){ //Compute double precision rsd and also new RHS vector. Linop_d.HermOp(sol_d, tmp_d); + if(shift) axpy(tmp_d,*shift,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; + ConjugateGradientShifted CG_d(Tolerance, MaxInnerIterations); + CG_d(Linop_d, src_d_in, sol_d,shift); TotalTimer.Stop(); - std::cout< &Linop, const Field &src, Field &psi) std::vector results(nshift,grid); (*this)(Linop,src,results,psi); } + void operator() (LinearOperatorBase &Linop, const Field &src, std::vector &results, Field &psi) { int nshift = shifts.order; @@ -105,11 +107,12 @@ void operator() (LinearOperatorBase &Linop, const Field &src, std::vector RealD a,b,c,d; RealD cp,bp,qq; //prev + int cb=src.checkerboard; // Matrix mult fields Field r(grid); - Field p(grid); + Field p(grid); p.checkerboard = src.checkerboard; Field tmp(grid); - Field mmp(grid); + Field mmp(grid);mmp.checkerboard = src.checkerboard; // Check lightest mass for(int s=0;s &Linop, const Field &src, std::vector p=src; //MdagM+m[0] + std::cout << "p.checkerboard " << p.checkerboard + << "mmp.checkerboard " << mmp.checkerboard << std::endl; + Linop.HermOpAndNorm(p,mmp,d,qq); axpy(mmp,mass[0],p,mmp); RealD rn = norm2(p); @@ -269,6 +275,7 @@ void operator() (LinearOperatorBase &Linop, const Field &src, std::vector RealD cn = norm2(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 */ +#ifndef GRID_CONJUGATE_GRADIENT_MULTI_MIXED_PREC_H +#define GRID_CONJUGATE_GRADIENT_MULTI_MIXED_PREC_H + +namespace Grid { + + //Mixed precision restarted defect correction CG + template::value == 2, int>::type = 0 +//, typename std::enable_if< getPrecision::value == 1, int>::type = 0 +> + class MixedPrecisionConjugateGradientMultiShift : public LinearFunction { + public: +// RealD Tolerance; + 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; + MultiShiftFunction shifts; + Integer iter; + + //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess +// LinearFunction *guesser; + + MixedPrecisionConjugateGradientMultiShift(GridBase* _sp_grid, LinearOperatorBase &_Linop_f, LinearOperatorBase &_Linop_d, +Integer maxinnerit, MultiShiftFunction &_shifts ) : + Linop_f(_Linop_f), Linop_d(_Linop_d), + MaxInnerIterations(maxinnerit), SinglePrecGrid(_sp_grid), + OuterLoopNormMult(100.), shifts(_shifts) {}; + + + void operator() (const FieldD &src_d_in, FieldD &sol_d){ + assert(0); // not yet implemented + } + void operator() (const FieldD &src_d_in, std::vector &sol_d){ + GridStopWatch TotalTimer; + TotalTimer.Start(); + + int cb = src_d_in.checkerboard; + + int nshift = shifts.order; + assert(nshift == sol_d.size()); + for(int i=0;i sol_f(nshift,SinglePrecGrid); + for(int i=0;i CG_f(inner_tol, MaxInnerIterations); + ConjugateGradientMultiShift MSCG(MaxInnerIterations,shifts); +// CG_f.ErrorOnNoConverge = false; + + GridStopWatch InnerCGtimer; + + GridStopWatch PrecChangeTimer; + +{ +// std::cout< &mass(shifts.poles); // Make references to array in "shifts" + std::vector &mresidual(shifts.tolerances); + std::vector alpha(nshift,1.); + std::vector ps(nshift,grid);// Search directions + + assert(sol_f.size()==nshift); + assert(mass.size()==nshift); + assert(mresidual.size()==nshift); + + // dynamic sized arrays on stack; 2d is a pain with vector + RealD bs[nshift]; + RealD rsq[nshift]; + RealD z[nshift][2]; + int converged[nshift]; + + const int primary =0; + + //Primary shift fields CG iteration + RealD a,b,c,d; + RealD cp,bp,qq; //prev + + int cb=src_f.checkerboard; + // Matrix mult fields + FieldF r(grid); r.checkerboard = src_f.checkerboard; + FieldF p(grid); p.checkerboard = src_f.checkerboard; + FieldF tmp(grid); tmp.checkerboard = src_f.checkerboard; + FieldF mmp(grid);mmp.checkerboard = src_f.checkerboard; + FieldF psi(grid);psi.checkerboard = src_f.checkerboard; + std::cout.precision(12); + std::cout< 2+Ls, so ~ 3x saving + // Pipelined CG gain: + // + // New Kernel: Load r, vector of coeffs, vector of pointers ps + // New Kernel: Load sol_f[0], vector of coeffs, vector of pointers ps + // If can predict the coefficient bs then we can fuse these and avoid write reread cyce + // on ps[s]. + // Before: 3 x npole + 3 x npole + // After : 2 x npole (ps[s]) => 3x speed up of multishift CG. + + if( (!converged[s]) ) { + axpy(sol_f[ss],-bs[s]*alpha[s],ps[s],sol_f[ss]); + } + } + + + if (k%MaxInnerIterations==0){ +// if (c < 1e-4*c_relup){ + RealD c_f=c; + precisionChange(tmp_d,psi); + RealD sol_norm =axpy_norm (psi_d,1.,tmp_d,psi_d); + tmp1 = norm2(psi); + zeroit(psi); + tmp2 = norm2(psi); + 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 */ +#ifndef GRID_CONJUGATE_GRADIENT_SHIFTED_H +#define GRID_CONJUGATE_GRADIENT_SHIFTED_H + +namespace Grid { + + ///////////////////////////////////////////////////////////// + // Base classes for iterative processes based on operators + // single input vec, single output vec. + ///////////////////////////////////////////////////////////// + + template + class ConjugateGradientShifted : public OperatorFunction { +public: + bool ErrorOnNoConverge; //throw an assert when the CG fails to converge. Defaults true. + RealD Tolerance; + Integer MaxIterations; + ConjugateGradientShifted(RealD tol,Integer maxit, bool err_on_no_conv = true) : Tolerance(tol), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv) { + }; + + void operator() (LinearOperatorBase &Linop,const Field &src, Field &psi ){ + (*this)(Linop,src,psi,NULL); + } + + void operator() (LinearOperatorBase &Linop,const Field &src, Field &psi, RealD *shift){ + + psi.checkerboard = src.checkerboard; + conformable(psi,src); + + RealD cp,c,a,d,b,ssq,qq,b_pred; + + Field p(src); + Field mmp(src); + Field r(src); + + //Initial residual computation & set up + RealD guess = norm2(psi); + assert(std::isnan(guess)==0); + + Linop.HermOpAndNorm(psi,mmp,d,b); + if(shift) axpy(mmp,*shift,psi,mmp); + RealD rn = norm2(psi); + if(shift) d += rn*(*shift); + RealD d2 = real(innerProduct(psi,mmp)); + b= norm2(mmp); + RealD src_norm=norm2(src); + r= src-mmp; + p= r; + + a =norm2(p); + cp =a; + ssq=norm2(src); + + std::cout< Author: paboyle +Author: Chulwoo Jung 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 diff --git a/lib/algorithms/iterative/Matrix.h b/lib/algorithms/iterative/Matrix.h index cb78a598..afdda2d5 100644 --- a/lib/algorithms/iterative/Matrix.h +++ b/lib/algorithms/iterative/Matrix.h @@ -36,7 +36,7 @@ Author: Peter Boyle #include #include #include -#include +#include /** Sign function **/