From 4180a4a8a710b3f699c0f16f32518b7dd296c421 Mon Sep 17 00:00:00 2001 From: David Murphy Date: Tue, 10 Dec 2019 17:20:35 -0500 Subject: [PATCH 1/2] Import BiCGSTAB solvers and tests --- Grid/algorithms/Algorithms.h | 2 + Grid/algorithms/LinearOperator.h | 132 +++++++++++ Grid/algorithms/iterative/BiCGSTAB.h | 221 ++++++++++++++++++ Grid/algorithms/iterative/BiCGSTABMixedPrec.h | 158 +++++++++++++ Grid/algorithms/iterative/SchurRedBlack.h | 135 +++++++++++ Grid/qcd/modules/Registration.h | 2 + Grid/qcd/modules/SolverModules.h | 11 + .../solver/Test_wilsonclover_bicgstab_prec.cc | 85 +++++++ .../Test_wilsonclover_bicgstab_schur.cc | 81 +++++++ .../Test_wilsonclover_bicgstab_unprec.cc | 80 +++++++ tests/solver/Test_wilsonclover_cg_prec.cc | 91 ++++++++ tests/solver/Test_wilsonclover_cg_schur.cc | 77 ++++++ tests/solver/Test_wilsonclover_cg_unprec.cc | 80 +++++++ .../Test_wilsonclover_mixedbicgstab_prec.cc | 126 ++++++++++ .../solver/Test_wilsonclover_mixedcg_prec.cc | 126 ++++++++++ 15 files changed, 1407 insertions(+) create mode 100644 Grid/algorithms/iterative/BiCGSTAB.h create mode 100644 Grid/algorithms/iterative/BiCGSTABMixedPrec.h create mode 100644 tests/solver/Test_wilsonclover_bicgstab_prec.cc create mode 100644 tests/solver/Test_wilsonclover_bicgstab_schur.cc create mode 100644 tests/solver/Test_wilsonclover_bicgstab_unprec.cc create mode 100644 tests/solver/Test_wilsonclover_cg_prec.cc create mode 100644 tests/solver/Test_wilsonclover_cg_schur.cc create mode 100644 tests/solver/Test_wilsonclover_cg_unprec.cc create mode 100644 tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc create mode 100644 tests/solver/Test_wilsonclover_mixedcg_prec.cc diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index f1ac1c81..b42338d5 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -41,11 +41,13 @@ Author: Peter Boyle #include #include +#include #include #include #include #include #include +#include #include #include #include diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index a768f431..4a5eeeec 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -320,6 +320,138 @@ public: return axpy_norm(out,-1.0,tmp,in); } }; + + template + class HermitianSchurOperatorBase : public LinearOperatorBase + { + public: + virtual RealD Mpc (const Field& in, Field& out) = 0; + virtual RealD MpcDag (const Field& in, Field& out) = 0; + virtual void MpcDagMpc(const Field& in, Field& out, RealD& ni, RealD& no) { + Field tmp(in.Grid()); + tmp.Checkerboard() = in.Checkerboard(); + ni = Mpc(in,tmp); + no = MpcDag(tmp,out); + } + virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) { + out.Checkerboard() = in.Checkerboard(); + + Mpc(in, out); + + ComplexD dot = innerProduct(in,out); n1 = real(dot); + n2 = norm2(out); + } + virtual void HermOp(const Field& in, Field& out) { + RealD n1, n2; + HermOpAndNorm(in, out, n1, n2); + } + void Op(const Field& in, Field& out) { + Mpc(in, out); + } + void AdjOp(const Field& in, Field& out) { + MpcDag(in, out); + } + // Support for coarsening to a multigrid + void OpDiag(const Field& in, Field& out) { + assert(0); // must coarsen the unpreconditioned system + } + void OpDir(const Field& in, Field& out, int dir, int disp) { + assert(0); + } + }; + + template + class HermitianSchurDiagMooeeOperator : public HermitianSchurOperatorBase + { + public: + Matrix& _Mat; + HermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; + virtual RealD Mpc(const Field& in, Field& out) { + Field tmp(in.Grid()); + tmp.Checkerboard() = !in.Checkerboard(); + + _Mat.Meooe(in, tmp); + _Mat.MooeeInv(tmp, out); + _Mat.Meooe(out, tmp); + + _Mat.Mooee(in, out); + + return axpy_norm(out, -1.0, tmp, out); + } + virtual RealD MpcDag(const Field& in, Field& out) { + Field tmp(in.Grid()); + + _Mat.MeooeDag(in, tmp); + _Mat.MooeeInvDag(tmp, out); + _Mat.MeooeDag(out, tmp); + + _Mat.MooeeDag(in, out); + + return axpy_norm(out, -1.0, tmp, out); + } + }; + + template + class HermitianSchurDiagOneOperator : public HermitianSchurOperatorBase + { + protected: + Matrix &_Mat; + + public: + HermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; + virtual RealD Mpc(const Field& in, Field& out) { + Field tmp(in.Grid()); + + _Mat.Meooe(in, out); + _Mat.MooeeInv(out, tmp); + _Mat.Meooe(tmp, out); + _Mat.MooeeInv(out, tmp); + + return axpy_norm(out, -1.0, tmp, in); + } + virtual RealD MpcDag(const Field& in, Field& out) { + Field tmp(in.Grid()); + + _Mat.MooeeInvDag(in, out); + _Mat.MeooeDag(out, tmp); + _Mat.MooeeInvDag(tmp, out); + _Mat.MeooeDag(out, tmp); + + return axpy_norm(out, -1.0, tmp, in); + } + }; + + template + class HermitianSchurDiagTwoOperator : public HermitianSchurOperatorBase + { + protected: + Matrix& _Mat; + + public: + HermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; + + virtual RealD Mpc(const Field& in, Field& out) { + Field tmp(in.Grid()); + + _Mat.MooeeInv(in, out); + _Mat.Meooe(out, tmp); + _Mat.MooeeInv(tmp, out); + _Mat.Meooe(out, tmp); + + return axpy_norm(out, -1.0, tmp, in); + } + virtual RealD MpcDag(const Field& in, Field& out) { + Field tmp(in.Grid()); + + _Mat.MeooeDag(in, out); + _Mat.MooeeInvDag(out, tmp); + _Mat.MeooeDag(tmp, out); + _Mat.MooeeInvDag(out, tmp); + + return axpy_norm(out, -1.0, tmp, in); + } + }; + /////////////////////////////////////////////////////////////////////////////////////////////////// // Left handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta --> ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta // Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta --> ( 1 - Moe Mee^-1 Meo ) Moo^-1 phi=eta ; psi = Moo^-1 phi diff --git a/Grid/algorithms/iterative/BiCGSTAB.h b/Grid/algorithms/iterative/BiCGSTAB.h new file mode 100644 index 00000000..3e609241 --- /dev/null +++ b/Grid/algorithms/iterative/BiCGSTAB.h @@ -0,0 +1,221 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/BiCGSTAB.h + +Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle +Author: paboyle +Author: juettner +Author: David Murphy + +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_BICGSTAB_H +#define GRID_BICGSTAB_H + +NAMESPACE_BEGIN(Grid); + +///////////////////////////////////////////////////////////// +// Base classes for iterative processes based on operators +// single input vec, single output vec. +///////////////////////////////////////////////////////////// + +template +class BiCGSTAB : public OperatorFunction +{ + public: + using OperatorFunction::operator(); + + bool ErrorOnNoConverge; // throw an assert when the CG fails to converge. + // Defaults true. + RealD Tolerance; + Integer MaxIterations; + Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion + + BiCGSTAB(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) + { + psi.Checkerboard() = src.Checkerboard(); + conformable(psi, src); + + RealD cp(0), rho(1), rho_prev(0), alpha(1), beta(0), omega(1); + RealD a(0), bo(0), d(0), b(0), ssq(0), qq(0); + + Field p(src); + Field r(src); + Field rhat(src); + Field v(src); + Field s(src); + Field t(src); + Field h(src); + + v = Zero(); + p = Zero(); + + // Initial residual computation & set up + RealD guess = norm2(psi); + assert(std::isnan(guess) == 0); + + Linop.HermOpAndNorm(psi, v, d, b); + + r = src - v; + rhat = r; + a = norm2(r); + ssq = norm2(src); + + std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: guess " << guess << std::endl; + std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: src " << ssq << std::endl; + std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: mp " << b << std::endl; + std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: r " << a << std::endl; + + RealD rsq = Tolerance * Tolerance * ssq; + + // Check if guess is really REALLY good :) + if(a <= rsq){ return; } + + std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: k=0 residual " << a << " target " << rsq << std::endl; + + GridStopWatch LinalgTimer; + GridStopWatch InnerTimer; + GridStopWatch AxpyNormTimer; + GridStopWatch LinearCombTimer; + GridStopWatch MatrixTimer; + GridStopWatch SolverTimer; + + SolverTimer.Start(); + int k; + for (k = 1; k <= MaxIterations; k++) + { + rho_prev = rho; + + LinalgTimer.Start(); + InnerTimer.Start(); + ComplexD Crho = innerProduct(rhat,r); + InnerTimer.Stop(); + rho = Crho.real(); + + beta = (rho / rho_prev) * (alpha / omega); + + LinearCombTimer.Start(); + bo = beta * omega; + auto p_v = p.View(); + auto r_v = r.View(); + auto v_v = v.View(); + accelerator_for(ss, p_v.size(), Field::vector_object::Nsimd(),{ + coalescedWrite(p_v[ss], beta*p_v(ss) - bo*v_v(ss) + r_v(ss)); + }); + LinearCombTimer.Stop(); + LinalgTimer.Stop(); + + MatrixTimer.Start(); + Linop.HermOp(p,v); + MatrixTimer.Stop(); + + LinalgTimer.Start(); + InnerTimer.Start(); + ComplexD Calpha = innerProduct(rhat,v); + InnerTimer.Stop(); + alpha = rho / Calpha.real(); + + LinearCombTimer.Start(); + auto h_v = h.View(); + auto psi_v = psi.View(); + accelerator_for(ss, h_v.size(), Field::vector_object::Nsimd(),{ + coalescedWrite(h_v[ss], alpha*p_v(ss) + psi_v(ss)); + }); + + auto s_v = s.View(); + accelerator_for(ss, s_v.size(), Field::vector_object::Nsimd(),{ + coalescedWrite(s_v[ss], -alpha*v_v(ss) + r_v(ss)); + }); + LinearCombTimer.Stop(); + LinalgTimer.Stop(); + + MatrixTimer.Start(); + Linop.HermOp(s,t); + MatrixTimer.Stop(); + + LinalgTimer.Start(); + InnerTimer.Start(); + ComplexD Comega = innerProduct(t,s); + InnerTimer.Stop(); + omega = Comega.real() / norm2(t); + + LinearCombTimer.Start(); + auto t_v = t.View(); + accelerator_for(ss, psi_v.size(), Field::vector_object::Nsimd(),{ + coalescedWrite(psi_v[ss], h_v(ss) + omega * s_v(ss)); + coalescedWrite(r_v[ss], -omega * t_v(ss) + s_v(ss)); + }); + LinearCombTimer.Stop(); + + cp = norm2(r); + LinalgTimer.Stop(); + + std::cout << GridLogIterative << "BiCGSTAB: Iteration " << k << " residual " << sqrt(cp/ssq) << " target " << Tolerance << std::endl; + + // Stopping condition + if(cp <= rsq) + { + SolverTimer.Stop(); + Linop.HermOpAndNorm(psi, v, d, qq); + p = v - src; + + RealD srcnorm = sqrt(norm2(src)); + RealD resnorm = sqrt(norm2(p)); + RealD true_residual = resnorm / srcnorm; + + std::cout << GridLogMessage << "BiCGSTAB Converged on iteration " << k << std::endl; + std::cout << GridLogMessage << "\tComputed residual " << sqrt(cp/ssq) << std::endl; + std::cout << GridLogMessage << "\tTrue residual " << true_residual << std::endl; + std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl; + + std::cout << GridLogMessage << "Time breakdown " << std::endl; + std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() << std::endl; + std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() << std::endl; + std::cout << GridLogMessage << "\tLinalg " << LinalgTimer.Elapsed() << std::endl; + std::cout << GridLogMessage << "\tInner " << InnerTimer.Elapsed() << std::endl; + std::cout << GridLogMessage << "\tAxpyNorm " << AxpyNormTimer.Elapsed() << std::endl; + std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() << std::endl; + + if(ErrorOnNoConverge){ assert(true_residual / Tolerance < 10000.0); } + + IterationsToComplete = k; + + return; + } + } + + std::cout << GridLogMessage << "BiCGSTAB did NOT converge" << std::endl; + + if(ErrorOnNoConverge){ assert(0); } + IterationsToComplete = k; + } +}; + +NAMESPACE_END(Grid); + +#endif diff --git a/Grid/algorithms/iterative/BiCGSTABMixedPrec.h b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h new file mode 100644 index 00000000..ebdb1071 --- /dev/null +++ b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h @@ -0,0 +1,158 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/algorithms/iterative/BiCGSTABMixedPrec.h + +Copyright (C) 2015 + +Author: Christopher Kelly +Author: David Murphy + +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_BICGSTAB_MIXED_PREC_H +#define GRID_BICGSTAB_MIXED_PREC_H + +NAMESPACE_BEGIN(Grid); + +// Mixed precision restarted defect correction BiCGSTAB +template::value == 2, int>::type = 0, typename std::enable_if< getPrecision::value == 1, int>::type = 0> +class MixedPrecisionBiCGSTAB : 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 + + //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess + LinearFunction *guesser; + + MixedPrecisionBiCGSTAB(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 operator() (const FieldD& src_d_in, FieldD& sol_d) + { + TotalInnerIterations = 0; + + GridStopWatch TotalTimer; + TotalTimer.Start(); + + int cb = src_d_in.Checkerboard(); + sol_d.Checkerboard() = cb; + + 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; + + 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 + + RealD inner_tol = InnerTolerance; + + FieldF src_f(SinglePrecGrid); + src_f.Checkerboard() = cb; + + FieldF sol_f(SinglePrecGrid); + sol_f.Checkerboard() = cb; + + BiCGSTAB 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++) + { + // 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 << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration " << outer_iter << " residual " << norm << " target " << stop << std::endl; + + if(norm < OuterLoopNormMult * stop){ + std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration converged on iteration " << outer_iter << std::endl; + break; + } + while(norm * inner_tol * inner_tol < stop){ inner_tol *= 2; } // inner_tol = sqrt(stop/norm) ?? + + PrecChangeTimer.Start(); + precisionChange(src_f, src_d); + PrecChangeTimer.Stop(); + + sol_f = Zero(); + + //Optionally improve inner solver guess (eg using known eigenvectors) + if(guesser != NULL){ (*guesser)(src_f, sol_f); } + + //Inner CG + CG_f.Tolerance = inner_tol; + InnerCGtimer.Start(); + CG_f(Linop_f, src_f, sol_f); + InnerCGtimer.Stop(); + TotalInnerIterations += CG_f.IterationsToComplete; + + //Convert sol back to double and add to double prec solution + PrecChangeTimer.Start(); + precisionChange(tmp_d, sol_f); + PrecChangeTimer.Stop(); + + axpy(sol_d, 1.0, tmp_d, sol_d); + } + + //Final trial CG + std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Starting final patch-up double-precision solve" << std::endl; + + BiCGSTAB CG_d(Tolerance, MaxInnerIterations); + CG_d(Linop_d, src_d_in, sol_d); + TotalFinalStepIterations = CG_d.IterationsToComplete; + + TotalTimer.Stop(); + std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl; + std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Total time " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl; + } +}; + +NAMESPACE_END(Grid); + +#endif diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index dd8a14b6..92f39cdd 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -405,6 +405,70 @@ namespace Grid { } }; + template class HermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase + { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + HermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + const bool _solnAsInitGuess = false) + : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + + ////////////////////////////////////////////////////// + // Override RedBlack specialisation + ////////////////////////////////////////////////////// + virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o) + { + GridBase* grid = _Matrix.RedBlackGrid(); + GridBase* fgrid = _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd , src_o, src); + + ///////////////////////////////////////////////////// + // src_o = Mdag * (source_o - Moe MeeInv source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInv(src_e, tmp); assert( tmp.Checkerboard() == Even ); + _Matrix.Meooe (tmp, Mtmp); assert( Mtmp.Checkerboard() == Odd ); + src_o -= Mtmp; assert( src_o.Checkerboard() == Odd ); + } + + virtual void RedBlackSolution(Matrix& _Matrix, const Field& sol_o, const Field& src_e, Field& sol) + { + GridBase* grid = _Matrix.RedBlackGrid(); + GridBase* fgrid = _Matrix.Grid(); + + Field tmp(grid); + Field sol_e(grid); + Field src_e_i(grid); + + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o, tmp); assert( tmp.Checkerboard() == Even ); + src_e_i = src_e - tmp; assert( src_e_i.Checkerboard() == Even ); + _Matrix.MooeeInv(src_e_i, sol_e); assert( sol_e.Checkerboard() == Even ); + + setCheckerboard(sol, sol_e); assert( sol_e.Checkerboard() == Even ); + setCheckerboard(sol, sol_o); assert( sol_o.Checkerboard() == Odd ); + } + + virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) + { + HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); assert(sol_o.Checkerboard() == Odd); + } + + virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) + { + HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + } + }; + /////////////////////////////////////////////////////////////////////////////////////////////////////// // Site diagonal is identity, right preconditioned by Mee^inv // ( 1 - Meo Moo^inv Moe Mee^inv ) phi =( 1 - Meo Moo^inv Moe Mee^inv ) Mee psi = = eta = eta @@ -482,5 +546,76 @@ namespace Grid { this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); } }; + + template class HermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase + { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + HermitianSchurRedBlackDiagTwoSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + const bool _solnAsInitGuess = false) + : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + + virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o) + { + GridBase* grid = _Matrix.RedBlackGrid(); + GridBase* fgrid = _Matrix.Grid(); + + Field tmp(grid); + Field Mtmp(grid); + + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd , src_o, src); + + ///////////////////////////////////////////////////// + // src_o = Mdag * (source_o - Moe MeeInv source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInv(src_e, tmp); assert( tmp.Checkerboard() == Even ); + _Matrix.Meooe (tmp, Mtmp); assert( Mtmp.Checkerboard() == Odd ); + src_o -= Mtmp; assert( src_o.Checkerboard() == Odd ); + } + + virtual void RedBlackSolution(Matrix& _Matrix, const Field& sol_o, const Field& src_e, Field& sol) + { + GridBase* grid = _Matrix.RedBlackGrid(); + GridBase* fgrid = _Matrix.Grid(); + + Field sol_o_i(grid); + Field tmp(grid); + Field sol_e(grid); + + //////////////////////////////////////////////// + // MooeeInv due to pecond + //////////////////////////////////////////////// + _Matrix.MooeeInv(sol_o, tmp); + sol_o_i = tmp; + + /////////////////////////////////////////////////// + // sol_e = M_ee^-1 * ( src_e - Meo sol_o )... + /////////////////////////////////////////////////// + _Matrix.Meooe(sol_o_i, tmp); assert( tmp.Checkerboard() == Even ); + tmp = src_e - tmp; assert( src_e.Checkerboard() == Even ); + _Matrix.MooeeInv(tmp, sol_e); assert( sol_e.Checkerboard() == Even ); + + setCheckerboard(sol, sol_e); assert( sol_e.Checkerboard() == Even ); + setCheckerboard(sol, sol_o_i); assert( sol_o_i.Checkerboard() == Odd ); + }; + + virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) + { + HermitianSchurDiagTwoOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + }; + + virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) + { + SchurDiagTwoOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + } + }; } + #endif diff --git a/Grid/qcd/modules/Registration.h b/Grid/qcd/modules/Registration.h index ebf8f420..ec28f020 100644 --- a/Grid/qcd/modules/Registration.h +++ b/Grid/qcd/modules/Registration.h @@ -80,6 +80,8 @@ static Registrar, static Registrar< ConjugateGradientModule, HMC_SolverModuleFactory > __CGWFmodXMLInit("ConjugateGradient"); +static Registrar< BiCGSTABModule, + HMC_SolverModuleFactory > __CGWFmodXMLInit("BiCGSTAB"); static Registrar< ConjugateResidualModule, HMC_SolverModuleFactory > __CRWFmodXMLInit("ConjugateResidual"); diff --git a/Grid/qcd/modules/SolverModules.h b/Grid/qcd/modules/SolverModules.h index 65cb91f9..2daaaed5 100644 --- a/Grid/qcd/modules/SolverModules.h +++ b/Grid/qcd/modules/SolverModules.h @@ -119,6 +119,17 @@ class ConjugateGradientModule: public SolverModule +class BiCGSTABModule: public SolverModule { + typedef SolverModule SolverBase; + using SolverBase::SolverBase; // for constructors + + // acquire resource + virtual void initialize(){ + this->SolverPtr.reset(new BiCGSTAB(this->Par_.tolerance, this->Par_.max_iterations, true)); + } +}; + template class ConjugateResidualModule: public SolverModule { typedef SolverModule SolverBase; diff --git a/tests/solver/Test_wilsonclover_bicgstab_prec.cc b/tests/solver/Test_wilsonclover_bicgstab_prec.cc new file mode 100644 index 00000000..708cc6f6 --- /dev/null +++ b/tests/solver/Test_wilsonclover_bicgstab_prec.cc @@ -0,0 +1,85 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_unprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=Zero(); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + double volume=1; + for(int mu=0;mu HermOp(Dw); + BiCGSTAB CG(1.0e-8,10000); + CG(HermOp, src_o, result_o); + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_bicgstab_schur.cc b/tests/solver/Test_wilsonclover_bicgstab_schur.cc new file mode 100644 index 00000000..c949748f --- /dev/null +++ b/tests/solver/Test_wilsonclover_bicgstab_schur.cc @@ -0,0 +1,81 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_unprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=Zero(); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + double volume=1; + for(int mu=0;mu CG(1.0e-8,10000); + HermitianSchurRedBlackDiagMooeeSolve SchurSolver(CG); + + SchurSolver(Dw, src, result); + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_bicgstab_unprec.cc b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc new file mode 100644 index 00000000..021b3247 --- /dev/null +++ b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc @@ -0,0 +1,80 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_unprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=Zero(); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + double volume=1; + for(int mu=0;mu HermOp(Dw); + BiCGSTAB CG(1.0e-8,10000); + CG(HermOp,src,result); + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_cg_prec.cc b/tests/solver/Test_wilsonclover_cg_prec.cc new file mode 100644 index 00000000..04202860 --- /dev/null +++ b/tests/solver/Test_wilsonclover_cg_prec.cc @@ -0,0 +1,91 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_prec.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=Zero(); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + std::vector U(4,&Grid); + + for(int mu=0;mu(Umu,mu); + } + + RealD mass = -0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + + // HermitianOperator HermOp(Dw); + // ConjugateGradient CG(1.0e-8,10000); + // CG(HermOp,src,result); + + LatticeFermion src_o(&RBGrid); + LatticeFermion result_o(&RBGrid); + pickCheckerboard(Odd,src_o,src); + result_o=Zero(); + + SchurDiagMooeeOperator HermOpEO(Dw); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOpEO,src_o,result_o); + + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_cg_schur.cc b/tests/solver/Test_wilsonclover_cg_schur.cc new file mode 100644 index 00000000..eaae24b3 --- /dev/null +++ b/tests/solver/Test_wilsonclover_cg_schur.cc @@ -0,0 +1,77 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_schur.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + LatticeFermion src(&Grid); random(pRNG,src); + LatticeFermion result(&Grid); result=Zero(); + LatticeFermion resid(&Grid); + + RealD mass = -0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + + ConjugateGradient CG(1.0e-8,10000); + SchurRedBlackDiagMooeeSolve SchurSolver(CG); + + SchurSolver(Dw,src,result); + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_cg_unprec.cc b/tests/solver/Test_wilsonclover_cg_unprec.cc new file mode 100644 index 00000000..49c52cdf --- /dev/null +++ b/tests/solver/Test_wilsonclover_cg_unprec.cc @@ -0,0 +1,80 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_wilson_cg_unprec.cc + + Copyright (C) 2015 + +Author: Azusa Yamaguchi +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + std::vector seeds({1,2,3,4}); + GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + + LatticeFermion src(&Grid); random(pRNG,src); + RealD nrm = norm2(src); + LatticeFermion result(&Grid); result=Zero(); + LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu); + + double volume=1; + for(int mu=0;mu HermOp(Dw); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOp,src,result); + + Grid_finalize(); +} diff --git a/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc new file mode 100644 index 00000000..848191c5 --- /dev/null +++ b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc @@ -0,0 +1,126 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_cg_prec.cc + + Copyright (C) 2015 + +Author: Peter Boyle + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::cout << GridLogMessage << "::::: NB: to enable a quick bit reproducibility check use the --checksums flag. " << std::endl; + + GridCartesian *FGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridCartesian *FGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d); + GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f); + + std::vector fSeeds({1, 2, 3, 4}); + GridParallelRNG fPRNG(FGrid_d); + fPRNG.SeedFixedIntegers(fSeeds); + + // clang-format off + LatticeFermionD src(FGrid_d); gaussian(fPRNG, src); + LatticeFermionD result(FGrid_d); result = Zero(); + LatticeGaugeFieldD Umu_d(FGrid_d); SU3::HotConfiguration(fPRNG, Umu_d); + LatticeGaugeFieldF Umu_f(FGrid_f); precisionChange(Umu_f, Umu_d); + // clang-format on + + RealD mass = -0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + WilsonCloverFermionD Dw_d(Umu_d, *FGrid_d, *FrbGrid_d, mass, csw_r, csw_t); + WilsonCloverFermionF Dw_f(Umu_f, *FGrid_f, *FrbGrid_f, mass, csw_r, csw_t); + + LatticeFermionD src_o(FrbGrid_d); + LatticeFermionD result_o(FrbGrid_d); + LatticeFermionD result_o_2(FrbGrid_d); + pickCheckerboard(Odd, src_o, src); + result_o.Checkerboard() = Odd; + result_o = Zero(); + result_o_2.Checkerboard() = Odd; + result_o_2 = Zero(); + + HermitianSchurDiagMooeeOperator HermOpEO_d(Dw_d); + HermitianSchurDiagMooeeOperator HermOpEO_f(Dw_f); + + std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl; + MixedPrecisionBiCGSTAB mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO_d); + mCG(src_o, result_o); + + std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; + BiCGSTAB CG(1.0e-8, 10000); + CG(HermOpEO_d, src_o, result_o_2); + + LatticeFermionD diff_o(FrbGrid_d); + RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2); + + std::cout << GridLogMessage << "::::::::::::: Diff between mixed and regular CG: " << diff << std::endl; + + #ifdef HAVE_LIME + if( GridCmdOptionExists(argv,argv+argc,"--checksums") ) + { + std::string file1("./Propagator1"); + emptyUserRecord record; + uint32_t nersc_csum; + uint32_t scidac_csuma; + uint32_t scidac_csumb; + typedef SpinColourVectorD FermionD; + typedef vSpinColourVectorD vFermionD; + + BinarySimpleMunger munge; + std::string format = getFormatString(); + + BinaryIO::writeLatticeObject(result_o,file1,munge, 0, format, + nersc_csum,scidac_csuma,scidac_csumb); + + std::cout << GridLogMessage << " Mixed checksums "<(result_o_2,file1,munge, 0, format, + nersc_csum,scidac_csuma,scidac_csumb); + + std::cout << GridLogMessage << " CG checksums "< + + 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 + +using namespace std; +using namespace Grid; + ; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + std::cout << GridLogMessage << "::::: NB: to enable a quick bit reproducibility check use the --checksums flag. " << std::endl; + + GridCartesian *FGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridCartesian *FGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d); + GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f); + + std::vector fSeeds({1, 2, 3, 4}); + GridParallelRNG fPRNG(FGrid_d); + fPRNG.SeedFixedIntegers(fSeeds); + + // clang-format off + LatticeFermionD src(FGrid_d); gaussian(fPRNG, src); + LatticeFermionD result(FGrid_d); result = Zero(); + LatticeGaugeFieldD Umu_d(FGrid_d); SU3::HotConfiguration(fPRNG, Umu_d); + LatticeGaugeFieldF Umu_f(FGrid_f); precisionChange(Umu_f, Umu_d); + // clang-format on + + RealD mass = -0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + WilsonCloverFermionD Dw_d(Umu_d, *FGrid_d, *FrbGrid_d, mass, csw_r, csw_t); + WilsonCloverFermionF Dw_f(Umu_f, *FGrid_f, *FrbGrid_f, mass, csw_r, csw_t); + + LatticeFermionD src_o(FrbGrid_d); + LatticeFermionD result_o(FrbGrid_d); + LatticeFermionD result_o_2(FrbGrid_d); + pickCheckerboard(Odd, src_o, src); + result_o.Checkerboard() = Odd; + result_o = Zero(); + result_o_2.Checkerboard() = Odd; + result_o_2 = Zero(); + + SchurDiagMooeeOperator HermOpEO_d(Dw_d); + SchurDiagMooeeOperator HermOpEO_f(Dw_f); + + std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl; + MixedPrecisionConjugateGradient mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO_d); + mCG(src_o, result_o); + + std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; + ConjugateGradient CG(1.0e-8, 10000); + CG(HermOpEO_d, src_o, result_o_2); + + LatticeFermionD diff_o(FrbGrid_d); + RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2); + + std::cout << GridLogMessage << "::::::::::::: Diff between mixed and regular CG: " << diff << std::endl; + + #ifdef HAVE_LIME + if( GridCmdOptionExists(argv,argv+argc,"--checksums") ) + { + std::string file1("./Propagator1"); + emptyUserRecord record; + uint32_t nersc_csum; + uint32_t scidac_csuma; + uint32_t scidac_csumb; + typedef SpinColourVectorD FermionD; + typedef vSpinColourVectorD vFermionD; + + BinarySimpleMunger munge; + std::string format = getFormatString(); + + BinaryIO::writeLatticeObject(result_o,file1,munge, 0, format, + nersc_csum,scidac_csuma,scidac_csumb); + + std::cout << GridLogMessage << " Mixed checksums "<(result_o_2,file1,munge, 0, format, + nersc_csum,scidac_csuma,scidac_csumb); + + std::cout << GridLogMessage << " CG checksums "< Date: Wed, 11 Dec 2019 11:46:18 -0500 Subject: [PATCH 2/2] Fix naming conventions to be consistent with Peter --- Grid/algorithms/LinearOperator.h | 24 ++++++---------- Grid/algorithms/iterative/BiCGSTAB.h | 11 ++++---- Grid/algorithms/iterative/BiCGSTABMixedPrec.h | 2 +- Grid/algorithms/iterative/SchurRedBlack.h | 28 +++++++++---------- .../solver/Test_wilsonclover_bicgstab_prec.cc | 2 +- .../Test_wilsonclover_bicgstab_schur.cc | 2 +- .../Test_wilsonclover_bicgstab_unprec.cc | 2 +- .../Test_wilsonclover_mixedbicgstab_prec.cc | 8 +++--- 8 files changed, 37 insertions(+), 42 deletions(-) diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 4a5eeeec..b70e2d6e 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -322,7 +322,7 @@ public: }; template - class HermitianSchurOperatorBase : public LinearOperatorBase + class NonHermitianSchurOperatorBase : public LinearOperatorBase { public: virtual RealD Mpc (const Field& in, Field& out) = 0; @@ -334,16 +334,10 @@ public: no = MpcDag(tmp,out); } virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) { - out.Checkerboard() = in.Checkerboard(); - - Mpc(in, out); - - ComplexD dot = innerProduct(in,out); n1 = real(dot); - n2 = norm2(out); + assert(0); } virtual void HermOp(const Field& in, Field& out) { - RealD n1, n2; - HermOpAndNorm(in, out, n1, n2); + assert(0); } void Op(const Field& in, Field& out) { Mpc(in, out); @@ -361,11 +355,11 @@ public: }; template - class HermitianSchurDiagMooeeOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase { public: Matrix& _Mat; - HermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); tmp.Checkerboard() = !in.Checkerboard(); @@ -392,13 +386,13 @@ public: }; template - class HermitianSchurDiagOneOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase { protected: Matrix &_Mat; public: - HermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); @@ -422,13 +416,13 @@ public: }; template - class HermitianSchurDiagTwoOperator : public HermitianSchurOperatorBase + class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase { protected: Matrix& _Mat; public: - HermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; + NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){}; virtual RealD Mpc(const Field& in, Field& out) { Field tmp(in.Grid()); diff --git a/Grid/algorithms/iterative/BiCGSTAB.h b/Grid/algorithms/iterative/BiCGSTAB.h index 3e609241..3a7be1ef 100644 --- a/Grid/algorithms/iterative/BiCGSTAB.h +++ b/Grid/algorithms/iterative/BiCGSTAB.h @@ -62,7 +62,7 @@ class BiCGSTAB : public OperatorFunction conformable(psi, src); RealD cp(0), rho(1), rho_prev(0), alpha(1), beta(0), omega(1); - RealD a(0), bo(0), d(0), b(0), ssq(0), qq(0); + RealD a(0), bo(0), b(0), ssq(0); Field p(src); Field r(src); @@ -79,7 +79,8 @@ class BiCGSTAB : public OperatorFunction RealD guess = norm2(psi); assert(std::isnan(guess) == 0); - Linop.HermOpAndNorm(psi, v, d, b); + Linop.Op(psi, v); + b = norm2(v); r = src - v; rhat = r; @@ -131,7 +132,7 @@ class BiCGSTAB : public OperatorFunction LinalgTimer.Stop(); MatrixTimer.Start(); - Linop.HermOp(p,v); + Linop.Op(p,v); MatrixTimer.Stop(); LinalgTimer.Start(); @@ -155,7 +156,7 @@ class BiCGSTAB : public OperatorFunction LinalgTimer.Stop(); MatrixTimer.Start(); - Linop.HermOp(s,t); + Linop.Op(s,t); MatrixTimer.Stop(); LinalgTimer.Start(); @@ -181,7 +182,7 @@ class BiCGSTAB : public OperatorFunction if(cp <= rsq) { SolverTimer.Stop(); - Linop.HermOpAndNorm(psi, v, d, qq); + Linop.Op(psi, v); p = v - src; RealD srcnorm = sqrt(norm2(src)); diff --git a/Grid/algorithms/iterative/BiCGSTABMixedPrec.h b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h index ebdb1071..d9b2a9d7 100644 --- a/Grid/algorithms/iterative/BiCGSTABMixedPrec.h +++ b/Grid/algorithms/iterative/BiCGSTABMixedPrec.h @@ -105,7 +105,7 @@ class MixedPrecisionBiCGSTAB : public LinearFunction 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); + Linop_d.Op(sol_d, tmp_d); RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration " << outer_iter << " residual " << norm << " target " << stop << std::endl; diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index 92f39cdd..d0b133a3 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -405,14 +405,14 @@ namespace Grid { } }; - template class HermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase + template class NonHermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase { public: typedef CheckerBoardedSparseMatrixBase Matrix; - HermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + NonHermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, const bool _solnAsInitGuess = false) - : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + : SchurRedBlackBase(RBSolver, initSubGuess, _solnAsInitGuess) {}; ////////////////////////////////////////////////////// // Override RedBlack specialisation @@ -458,14 +458,14 @@ namespace Grid { virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) { - HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); assert(sol_o.Checkerboard() == Odd); + NonHermitianSchurDiagMooeeOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); assert(sol_o.Checkerboard() == Odd); } virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) { - HermitianSchurDiagMooeeOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagMooeeOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); } }; @@ -547,7 +547,7 @@ namespace Grid { } }; - template class HermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase + template class NonHermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase { public: typedef CheckerBoardedSparseMatrixBase Matrix; @@ -555,9 +555,9 @@ namespace Grid { ///////////////////////////////////////////////////// // Wrap the usual normal equations Schur trick ///////////////////////////////////////////////////// - HermitianSchurRedBlackDiagTwoSolve(OperatorFunction& HermitianRBSolver, const bool initSubGuess = false, + NonHermitianSchurRedBlackDiagTwoSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, const bool _solnAsInitGuess = false) - : SchurRedBlackBase(HermitianRBSolver, initSubGuess, _solnAsInitGuess) {}; + : SchurRedBlackBase(RBSolver, initSubGuess, _solnAsInitGuess) {}; virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o) { @@ -606,14 +606,14 @@ namespace Grid { virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o) { - HermitianSchurDiagTwoOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); }; virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) { - SchurDiagTwoOperator _HermOpEO(_Matrix); - this->_HermitianRBSolver(_HermOpEO, src_o, sol_o); + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); } }; } diff --git a/tests/solver/Test_wilsonclover_bicgstab_prec.cc b/tests/solver/Test_wilsonclover_bicgstab_prec.cc index 708cc6f6..c1905400 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_prec.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_prec.cc @@ -77,7 +77,7 @@ int main (int argc, char ** argv) pickCheckerboard(Odd, src_o, src); result_o = Zero(); - HermitianSchurDiagMooeeOperator HermOp(Dw); + NonHermitianSchurDiagMooeeOperator HermOp(Dw); BiCGSTAB CG(1.0e-8,10000); CG(HermOp, src_o, result_o); diff --git a/tests/solver/Test_wilsonclover_bicgstab_schur.cc b/tests/solver/Test_wilsonclover_bicgstab_schur.cc index c949748f..af021267 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_schur.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_schur.cc @@ -73,7 +73,7 @@ int main (int argc, char ** argv) WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); BiCGSTAB CG(1.0e-8,10000); - HermitianSchurRedBlackDiagMooeeSolve SchurSolver(CG); + NonHermitianSchurRedBlackDiagMooeeSolve SchurSolver(CG); SchurSolver(Dw, src, result); diff --git a/tests/solver/Test_wilsonclover_bicgstab_unprec.cc b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc index 021b3247..21b7521d 100644 --- a/tests/solver/Test_wilsonclover_bicgstab_unprec.cc +++ b/tests/solver/Test_wilsonclover_bicgstab_unprec.cc @@ -72,7 +72,7 @@ int main (int argc, char ** argv) RealD csw_t = 1.0; WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); - HermitianLinearOperator HermOp(Dw); + NonHermitianLinearOperator HermOp(Dw); BiCGSTAB CG(1.0e-8,10000); CG(HermOp,src,result); diff --git a/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc index 848191c5..0af83f8b 100644 --- a/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc +++ b/tests/solver/Test_wilsonclover_mixedbicgstab_prec.cc @@ -80,16 +80,16 @@ int main (int argc, char ** argv) result_o_2.Checkerboard() = Odd; result_o_2 = Zero(); - HermitianSchurDiagMooeeOperator HermOpEO_d(Dw_d); - HermitianSchurDiagMooeeOperator HermOpEO_f(Dw_f); + NonHermitianSchurDiagMooeeOperator NonHermOpEO_d(Dw_d); + NonHermitianSchurDiagMooeeOperator NonHermOpEO_f(Dw_f); std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl; - MixedPrecisionBiCGSTAB mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO_d); + MixedPrecisionBiCGSTAB mCG(1.0e-8, 10000, 50, FrbGrid_f, NonHermOpEO_f, NonHermOpEO_d); mCG(src_o, result_o); std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; BiCGSTAB CG(1.0e-8, 10000); - CG(HermOpEO_d, src_o, result_o_2); + CG(NonHermOpEO_d, src_o, result_o_2); LatticeFermionD diff_o(FrbGrid_d); RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);