diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 97ab4dc1..f73a2398 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -42,11 +42,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 402d1244..50600d2d 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -334,6 +334,132 @@ public: return axpy_norm(out,-1.0,tmp,in); } }; + + template + class NonHermitianSchurOperatorBase : 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) { + assert(0); + } + virtual void HermOp(const Field& in, Field& out) { + assert(0); + } + 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 NonHermitianSchurDiagMooeeOperator : public NonHermitianSchurOperatorBase + { + public: + Matrix& _Mat; + NonHermitianSchurDiagMooeeOperator(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 NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase + { + protected: + Matrix &_Mat; + + public: + NonHermitianSchurDiagOneOperator (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 NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase + { + protected: + Matrix& _Mat; + + public: + NonHermitianSchurDiagTwoOperator(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..3a7be1ef --- /dev/null +++ b/Grid/algorithms/iterative/BiCGSTAB.h @@ -0,0 +1,222 @@ +/************************************************************************************* + +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), b(0), ssq(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.Op(psi, v); + b = norm2(v); + + 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.Op(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.Op(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.Op(psi, v); + 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..d9b2a9d7 --- /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.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; + + 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..d0b133a3 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -405,6 +405,70 @@ namespace Grid { } }; + template class NonHermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase + { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + NonHermitianSchurRedBlackDiagMooeeSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, + const bool _solnAsInitGuess = false) + : SchurRedBlackBase(RBSolver, 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) + { + 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) + { + NonHermitianSchurDiagMooeeOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, 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 NonHermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase + { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + ///////////////////////////////////////////////////// + // Wrap the usual normal equations Schur trick + ///////////////////////////////////////////////////// + NonHermitianSchurRedBlackDiagTwoSolve(OperatorFunction& RBSolver, const bool initSubGuess = false, + const bool _solnAsInitGuess = false) + : SchurRedBlackBase(RBSolver, 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) + { + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, src_o, sol_o); + }; + + virtual void RedBlackSolve(Matrix& _Matrix, const std::vector& src_o, std::vector& sol_o) + { + NonHermitianSchurDiagTwoOperator _OpEO(_Matrix); + this->_HermitianRBSolver(_OpEO, 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..c1905400 --- /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..af021267 --- /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); + NonHermitianSchurRedBlackDiagMooeeSolve 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..21b7521d --- /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..0af83f8b --- /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(); + + 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, NonHermOpEO_f, NonHermOpEO_d); + mCG(src_o, result_o); + + std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl; + BiCGSTAB CG(1.0e-8, 10000); + 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); + + 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 "<