From 48edb8f72ed73c8ed5f44ad9f13bbb656b6e97da Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 31 Mar 2021 19:31:46 +0200 Subject: [PATCH 01/24] HMC prep for DDHMC --- Grid/qcd/hmc/HMC.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index 44674ea5..b05e0561 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -34,6 +34,7 @@ directory * @brief Classes for Hybrid Monte Carlo update * * @author Guido Cossu + * @author Peter Boyle */ //-------------------------------------------------------------------- #pragma once From aaf5ebf34540920f8a1e16a32fd648defbe20636 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 31 Mar 2021 19:32:57 +0200 Subject: [PATCH 02/24] Small hack --- HMC/Mobius2p1fRHMC.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/HMC/Mobius2p1fRHMC.cc b/HMC/Mobius2p1fRHMC.cc index b958d548..91bf7601 100644 --- a/HMC/Mobius2p1fRHMC.cc +++ b/HMC/Mobius2p1fRHMC.cc @@ -56,12 +56,12 @@ int main(int argc, char **argv) { MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 0; + HMCparams.StartTrajectory = 49; HMCparams.Trajectories = 200; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; - HMCparams.StartingType =std::string("ColdStart"); - // HMCparams.StartingType =std::string("CheckpointStart"); + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.MD = MD; HMCWrapper TheHMC(HMCparams); @@ -71,7 +71,7 @@ int main(int argc, char **argv) { CheckpointerParameters CPparams; CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.rng_prefix = "ckpoint_EODWF_rng"; - CPparams.saveInterval = 10; + CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); From ccd30e14851a027e639839d67b1fb674296baffd Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:03:01 +0200 Subject: [PATCH 03/24] 4D pseudofermion in Cayley action --- .../TwoFlavourRatio4DPseudoFermion.h | 197 ++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h b/Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h new file mode 100644 index 00000000..f87acee9 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h @@ -0,0 +1,197 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + + +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////// +// Two flavour ratio +/////////////////////////////////////// +template +class TwoFlavourRatio4DPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField phi4; // the pseudo fermion field for this trajectory + +public: + TwoFlavourRatio4DPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS + ) : NumOp(_NumOp), + DenOp(_DenOp), + DerivativeSolver(DS), + ActionSolver(AS), + phi4(_NumOp.GaugeGrid()) + {}; + + virtual std::string action_name(){return "TwoFlavourRatio4DPseudoFermionAction";} + + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta4(NumOp.GaugeGrid()); + FermionField eta5(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); + FermionField phi5(NumOp.FermionGrid()); + + gaussian(pRNG,eta4); + NumOp.ImportFourDimPseudoFermion(eta4,eta5); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + + DenOp.M(eta5,phi5); // M eta + NumOp.Mdag(phi5,tmp); // Vdag M eta + phi5 = Zero(); + ActionSolver(MdagMOp,tmp,phi5); // (VdagV)^-1 M eta = V^-1 Vdag^-1 Vdag M eta = V^-1 M eta + phi5=phi5*scale; + + // Project to 4d + NumOp.ExportFourDimPseudoFermion(phi5,phi4); + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField Y4(NumOp.GaugeGrid()); + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField phi5(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + NumOp.ImportFourDimPseudoFermion(phi4,phi5); + NumOp.M(phi5,Y); // Y= V phi + DenOp.Mdag(Y,X); // X= Mdag V phi + Y=Zero(); + ActionSolver(MdagMOp,X,Y); // Y= (MdagM)^-1 Mdag Vdag phi = M^-1 V phi + + NumOp.ExportFourDimPseudoFermion(Y,Y4); + + RealD action = norm2(Y4); + + return action; + }; + + ////////////////////////////////////////////////////// + // dS/du = 2 Re phi^dag (V^dag M^-dag)_11 (M^-1 d V)_11 phi + // - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField phi(NumOp.FermionGrid()); + FermionField Vphi(NumOp.FermionGrid()); + FermionField MinvVphi(NumOp.FermionGrid()); + FermionField tmp4(NumOp.GaugeGrid()); + FermionField MdagInvMinvVphi(NumOp.FermionGrid()); + + GaugeField force(NumOp.GaugeGrid()); + + //Y=V phi + //X = (Mdag V phi + //Y = (Mdag M)^-1 Mdag V phi = M^-1 V Phi + NumOp.ImportFourDimPseudoFermion(phi4,phi); + NumOp.M(phi,Vphi); // V phi + DenOp.Mdag(Vphi,X); // X= Mdag V phi + Y=Zero(); + DerivativeSolver(MdagMOp,X,MinvVphi);// M^-1 V phi + + // Projects onto the physical space and back + NumOp.ExportFourDimPseudoFermion(MinvVphi,tmp4); + NumOp.ImportFourDimPseudoFermion(tmp4,Y); + + X=Zero(); + DerivativeSolver(MdagMOp,Y,X);// X = (MdagM)^-1 proj M^-1 V phi + DenOp.M(X,MdagInvMinvVphi); + + // phi^dag (Vdag Mdag^-1) (M^-1 dV) phi + NumOp.MDeriv(force ,MdagInvMinvVphi , phi, DaggerNo ); dSdU=force; + + // phi^dag (dVdag Mdag^-1) (M^-1 V) phi + NumOp.MDeriv(force , phi, MdagInvMinvVphi ,DaggerYes ); dSdU=dSdU+force; + + // - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi + DenOp.MDeriv(force,MdagInvMinvVphi,MinvVphi,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,MinvVphi,MdagInvMinvVphi,DaggerYes); dSdU=dSdU-force; + + dSdU *= -1.0; + //dSdU = - Ta(dSdU); + + }; +}; + +NAMESPACE_END(Grid); + + From a6326b664e1499cca7be33ee6e4741ccf99f7b6f Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:04:16 +0200 Subject: [PATCH 04/24] Move momenutm filter to earlier in the include sequence so it can be used by DDHMC --- Grid/qcd/action/ActionCore.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Grid/qcd/action/ActionCore.h b/Grid/qcd/action/ActionCore.h index 6544318d..2ec74a95 100644 --- a/Grid/qcd/action/ActionCore.h +++ b/Grid/qcd/action/ActionCore.h @@ -58,6 +58,8 @@ NAMESPACE_CHECK(Scalar); //////////////////////////////////////////// // Utility functions //////////////////////////////////////////// +#include +//#include #include NAMESPACE_CHECK(Metric); #include From 0765f303081bf4eaeb5d1018ac5b59ca48428524 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:05:42 +0200 Subject: [PATCH 05/24] 4D pseudo fermion support --- Grid/qcd/action/fermion/CayleyFermion5D.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index c7d68d73..54608137 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -60,6 +60,8 @@ public: /////////////////////////////////////////////////////////////// virtual void Dminus(const FermionField &psi, FermionField &chi); virtual void DminusDag(const FermionField &psi, FermionField &chi); + virtual void ImportFourDimPseudoFermion(const FermionField &input,FermionField &imported); + virtual void ExportFourDimPseudoFermion(const FermionField &solution,FermionField &exported); virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d); virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d); virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d); From fe00c964355842754b8bbb0a4c58c9cbccb79d09 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:06:11 +0200 Subject: [PATCH 06/24] 4D Pseudofermion support --- Grid/qcd/action/fermion/FermionOperator.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/Grid/qcd/action/fermion/FermionOperator.h b/Grid/qcd/action/fermion/FermionOperator.h index 570e350d..0f379d30 100644 --- a/Grid/qcd/action/fermion/FermionOperator.h +++ b/Grid/qcd/action/fermion/FermionOperator.h @@ -171,6 +171,16 @@ public: /////////////////////////////////////////////// virtual void Dminus(const FermionField &psi, FermionField &chi) { chi=psi; } virtual void DminusDag(const FermionField &psi, FermionField &chi) { chi=psi; } + + virtual void ImportFourDimPseudoFermion(const FermionField &input,FermionField &imported) + { + imported = input; + }; + virtual void ExportFourDimPseudoFermion(const FermionField &solution,FermionField &exported) + { + exported=solution; + }; + virtual void ImportPhysicalFermionSource(const FermionField &input,FermionField &imported) { imported = input; From 09288d633c9897916e20b85e70dccbf2b7a33733 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:06:52 +0200 Subject: [PATCH 07/24] 4D pseudofermoin --- .../CayleyFermion5DImplementation.h | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h index f11e9c44..8c01e064 100644 --- a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h @@ -112,7 +112,6 @@ void CayleyFermion5D::ImportUnphysicalFermion(const FermionField &input4d, axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); imported5d=tmp; } - template void CayleyFermion5D::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d) { @@ -127,6 +126,37 @@ void CayleyFermion5D::ImportPhysicalFermionSource(const FermionField &inpu axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); Dminus(tmp,imported5d); } +//////////////////////////////////////////////////// +// Added for fourD pseudofermion det estimation +//////////////////////////////////////////////////// +template +void CayleyFermion5D::ImportFourDimPseudoFermion(const FermionField &input4d,FermionField &imported5d) +{ + int Ls = this->Ls; + FermionField tmp(this->FermionGrid()); + conformable(imported5d.Grid(),this->FermionGrid()); + conformable(input4d.Grid() ,this->GaugeGrid()); + tmp = Zero(); + InsertSlice(input4d, tmp, 0 , 0); + InsertSlice(input4d, tmp, Ls-1, 0); + axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, 0, 0); + axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); + imported5d=tmp; +} +template +void CayleyFermion5D::ExportFourDimPseudoFermion(const FermionField &solution5d,FermionField &exported4d) +{ + int Ls = this->Ls; + FermionField tmp(this->FermionGrid()); + tmp = solution5d; + conformable(solution5d.Grid(),this->FermionGrid()); + conformable(exported4d.Grid(),this->GaugeGrid()); + axpby_ssp_pminus(tmp, 0., solution5d, 1., solution5d, 0, 0); + axpby_ssp_pplus (tmp, 1., tmp , 1., solution5d, 0, Ls-1); + ExtractSlice(exported4d, tmp, 0, 0); +} + +// Dminus template void CayleyFermion5D::Dminus(const FermionField &psi, FermionField &chi) { From 21165ed48961bca438136980ecceb86724be8ef1 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Sat, 10 Apr 2021 01:08:23 +0200 Subject: [PATCH 08/24] Better logging and moving momentumfilter --- Grid/qcd/hmc/HMC.h | 80 ++++++++++++++----- Grid/qcd/hmc/integrators/Integrator.h | 1 - Grid/qcd/hmc/integrators/MomentumFilter.h | 94 ----------------------- 3 files changed, 59 insertions(+), 116 deletions(-) delete mode 100644 Grid/qcd/hmc/integrators/MomentumFilter.h diff --git a/Grid/qcd/hmc/HMC.h b/Grid/qcd/hmc/HMC.h index b05e0561..2e391caa 100644 --- a/Grid/qcd/hmc/HMC.h +++ b/Grid/qcd/hmc/HMC.h @@ -116,22 +116,17 @@ private: random(sRNG, rn_test); - std::cout << GridLogMessage - << "--------------------------------------------------\n"; - std::cout << GridLogMessage << "exp(-dH) = " << prob - << " Random = " << rn_test << "\n"; - std::cout << GridLogMessage - << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "exp(-dH) = " << prob << " Random = " << rn_test << "\n"; + std::cout << GridLogMessage << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n"; if ((prob > 1.0) || (rn_test <= prob)) { // accepted std::cout << GridLogMessage << "Metropolis_test -- ACCEPTED\n"; - std::cout << GridLogMessage - << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; return true; } else { // rejected std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n"; - std::cout << GridLogMessage - << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; return false; } } @@ -140,18 +135,63 @@ private: // Evolution ///////////////////////////////////////////////////////// RealD evolve_hmc_step(Field &U) { - TheIntegrator.refresh(U, sRNG, pRNG); // set U and initialize P and phi's - RealD H0 = TheIntegrator.S(U); // initial state action + GridBase *Grid = U.Grid(); + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // Mainly for DDHMC perform a random translation of U modulo volume + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Random shifting gauge field by ["; + for(int d=0;dNd();d++) { + + int L = Grid->GlobalDimensions()[d]; + + RealD rn_uniform; random(sRNG, rn_uniform); + + int shift = (int) (rn_uniform*L); + + std::cout << shift; + if(dNd()-1) std::cout <<","; + else std::cout <<"]\n"; + + U = Cshift(U,d,shift); + } + std::cout << GridLogMessage << "--------------------------------------------------\n"; + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // set U and initialize P and phi's + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Refresh momenta and pseudofermions"; + TheIntegrator.refresh(U, sRNG, pRNG); + std::cout << GridLogMessage << "--------------------------------------------------\n"; + + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // initial state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Compute initial action"; + RealD H0 = TheIntegrator.S(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; std::streamsize current_precision = std::cout.precision(); std::cout.precision(15); std::cout << GridLogMessage << "Total H before trajectory = " << H0 << "\n"; std::cout.precision(current_precision); + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << " Molecular Dynamics evolution "; TheIntegrator.integrate(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; - RealD H1 = TheIntegrator.S(U); // updated state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + // updated state action + ////////////////////////////////////////////////////////////////////////////////////////////////////// + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Compute final action"; + RealD H1 = TheIntegrator.S(U); + std::cout << GridLogMessage << "--------------------------------------------------\n"; /////////////////////////////////////////////////////////// if(0){ @@ -164,17 +204,14 @@ private: } /////////////////////////////////////////////////////////// - std::cout.precision(15); - std::cout << GridLogMessage << "Total H after trajectory = " << H1 - << " dH = " << H1 - H0 << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; + std::cout << GridLogMessage << "Total H after trajectory = " << H1 << " dH = " << H1 - H0 << "\n"; + std::cout << GridLogMessage << "--------------------------------------------------\n"; std::cout.precision(current_precision); return (H1 - H0); } - - - public: ///////////////////////////////////////// @@ -196,8 +233,11 @@ public: // Actual updates (evolve a copy Ucopy then copy back eventually) unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory; + for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) { + std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n"; + if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) { std::cout << GridLogMessage << "-- Thermalization" << std::endl; } @@ -217,8 +257,6 @@ public: if (accept) Ucur = Ucopy; - - double t1=usecond(); std::cout << GridLogMessage << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl; diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index aa28c6c8..6b0e3caf 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -33,7 +33,6 @@ directory #define INTEGRATOR_INCLUDED #include -#include "MomentumFilter.h" NAMESPACE_BEGIN(Grid); diff --git a/Grid/qcd/hmc/integrators/MomentumFilter.h b/Grid/qcd/hmc/integrators/MomentumFilter.h deleted file mode 100644 index 2a15d80c..00000000 --- a/Grid/qcd/hmc/integrators/MomentumFilter.h +++ /dev/null @@ -1,94 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h - -Copyright (C) 2015 - -Author: Christopher Kelly -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 */ -//-------------------------------------------------------------------- -#ifndef MOMENTUM_FILTER -#define MOMENTUM_FILTER - -NAMESPACE_BEGIN(Grid); - -//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh - -template -struct MomentumFilterBase{ - virtual void applyFilter(MomentaField &P) const; -}; - -//Do nothing -template -struct MomentumFilterNone: public MomentumFilterBase{ - void applyFilter(MomentaField &P) const override{} -}; - -//Multiply each site/direction by a Lorentz vector complex number field -//Can be used to implement a mask, zeroing out sites -template -struct MomentumFilterApplyPhase: public MomentumFilterBase{ - typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type - typedef typename MomentaField::scalar_type scalar_type; //scalar complex type - typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction - typedef Lattice LatticeLorentzScalarType; - - LatticeLorentzScalarType phase; - - MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){} - - //Default to uniform field of (1,0) - MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){ - LorentzScalarType one; - for(int mu=0;mu Date: Thu, 6 May 2021 23:10:36 +0200 Subject: [PATCH 09/24] Mdagger solve support --- Grid/algorithms/LinearOperator.h | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/Grid/algorithms/LinearOperator.h b/Grid/algorithms/LinearOperator.h index 1add212c..5a73f881 100644 --- a/Grid/algorithms/LinearOperator.h +++ b/Grid/algorithms/LinearOperator.h @@ -223,9 +223,14 @@ class SchurOperatorBase : public LinearOperatorBase { Mpc(in,tmp); MpcDag(tmp,out); } + virtual void MpcMpcDag(const Field &in, Field &out) { + Field tmp(in.Grid()); + tmp.Checkerboard() = in.Checkerboard(); + MpcDag(in,tmp); + Mpc(tmp,out); + } virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ - out.Checkerboard() = in.Checkerboard(); - MpcDagMpc(in,out); + HermOp(in,out); ComplexD dot= innerProduct(in,out); n1=real(dot); n2=norm2(out); @@ -276,6 +281,16 @@ template axpy(out,-1.0,tmp,out); } }; +// Mpc MpcDag system presented as the HermOp +template +class SchurDiagMooeeDagOperator : public SchurDiagMooeeOperator { + public: + virtual void HermOp(const Field &in, Field &out){ + out.Checkerboard() = in.Checkerboard(); + this->MpcMpcDag(in,out); + } + SchurDiagMooeeDagOperator (Matrix &Mat): SchurDiagMooeeOperator(Mat){}; +}; template class SchurDiagOneOperator : public SchurOperatorBase { protected: From 45440da79ddde3c89e07bb2219ebdd82db2af702 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:12:34 +0200 Subject: [PATCH 10/24] Tuning integrator --- HMC/Mobius2p1fRHMC.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/HMC/Mobius2p1fRHMC.cc b/HMC/Mobius2p1fRHMC.cc index 91bf7601..00e14710 100644 --- a/HMC/Mobius2p1fRHMC.cc +++ b/HMC/Mobius2p1fRHMC.cc @@ -52,11 +52,11 @@ int main(int argc, char **argv) { // MD.name = std::string("Force Gradient"); typedef GenericHMCRunner HMCWrapper; MD.name = std::string("MinimumNorm2"); - MD.MDsteps = 20; + MD.MDsteps = 12; MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 49; + HMCparams.StartTrajectory = 139; HMCparams.Trajectories = 200; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; @@ -130,7 +130,7 @@ int main(int argc, char **argv) { // Collect actions //////////////////////////////////// ActionLevel Level1(1); - ActionLevel Level2(4); + ActionLevel Level2(8); //////////////////////////////////// // Strange action From cb285681989768016f86e91492230d75f487ba90 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:12:57 +0200 Subject: [PATCH 11/24] Tuning integrator --- HMC/Mobius2p1fEOFA.cc | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/HMC/Mobius2p1fEOFA.cc b/HMC/Mobius2p1fEOFA.cc index b1294da5..68a3a5b5 100644 --- a/HMC/Mobius2p1fEOFA.cc +++ b/HMC/Mobius2p1fEOFA.cc @@ -188,15 +188,15 @@ int main(int argc, char **argv) { IntegratorParameters MD; // typedef GenericHMCRunner HMCWrapper; // MD.name = std::string("Leap Frog"); - typedef GenericHMCRunner HMCWrapper; - MD.name = std::string("Force Gradient"); - // typedef GenericHMCRunner HMCWrapper; - // MD.name = std::string("MinimumNorm2"); - MD.MDsteps = 6; + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Force Gradient"); + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 12; MD.trajL = 1.0; HMCparameters HMCparams; - HMCparams.StartTrajectory = 590; + HMCparams.StartTrajectory = 211; HMCparams.Trajectories = 1000; HMCparams.NoMetropolisUntil= 0; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; @@ -209,9 +209,9 @@ int main(int argc, char **argv) { TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition CheckpointerParameters CPparams; - CPparams.config_prefix = "ckpoint_EODWF_lat"; - CPparams.rng_prefix = "ckpoint_EODWF_rng"; - CPparams.saveInterval = 10; + CPparams.config_prefix = "ckpoint_EOFA_lat"; + CPparams.rng_prefix = "ckpoint_EOFA_rng"; + CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); @@ -263,7 +263,7 @@ int main(int argc, char **argv) { FermionActionF::ImplParams ParamsF(boundary); double ActionStoppingCondition = 1e-10; - double DerivativeStoppingCondition = 1e-6; + double DerivativeStoppingCondition = 1e-8; double MaxCGIterations = 30000; //////////////////////////////////// From cdeb7182297b557d465950fcd3203a0828f94d53 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:15:16 +0200 Subject: [PATCH 12/24] 4D pseudo fermion, with Schur red black solvers --- .../TwoFlavourRatioEO4DPseudoFermion.h | 203 ++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 Grid/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h diff --git a/Grid/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h b/Grid/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h new file mode 100644 index 00000000..44e7b2d7 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h @@ -0,0 +1,203 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////// +// Two flavour ratio +/////////////////////////////////////// +template +class TwoFlavourRatioEO4DPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + typedef FermionOperator FermOp; + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &DerivativeDagSolver; + OperatorFunction &ActionSolver; + OperatorFunction &HeatbathSolver; + + FermionField phi4; // the pseudo fermion field for this trajectory + +public: + TwoFlavourRatioEO4DPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & AS ) : + TwoFlavourRatioEO4DPseudoFermionAction(_NumOp,_DenOp, DS,DS,AS,AS) {}; + TwoFlavourRatioEO4DPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + OperatorFunction & DS, + OperatorFunction & DDS, + OperatorFunction & AS, + OperatorFunction & HS + ) : NumOp(_NumOp), + DenOp(_DenOp), + DerivativeSolver(DS), + DerivativeDagSolver(DDS), + ActionSolver(AS), + HeatbathSolver(HS), + phi4(_NumOp.GaugeGrid()) + {}; + + virtual std::string action_name(){return "TwoFlavourRatioEO4DPseudoFermionAction";} + + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta4(NumOp.GaugeGrid()); + FermionField eta5(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); + FermionField phi5(NumOp.FermionGrid()); + + gaussian(pRNG,eta4); + NumOp.ImportFourDimPseudoFermion(eta4,eta5); + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + SchurRedBlackDiagMooeeSolve PrecSolve(HeatbathSolver); + + DenOp.M(eta5,tmp); // M eta + PrecSolve(NumOp,tmp,phi5); // phi = V^-1 M eta + phi5=phi5*scale; + std::cout << GridLogMessage << "4d pf refresh "<< norm2(phi5)<<"\n"; + // Project to 4d + NumOp.ExportFourDimPseudoFermion(phi5,phi4); + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField Y4(NumOp.GaugeGrid()); + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField phi5(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + SchurRedBlackDiagMooeeSolve PrecSolve(ActionSolver); + + NumOp.ImportFourDimPseudoFermion(phi4,phi5); + NumOp.M(phi5,X); // X= V phi + PrecSolve(DenOp,X,Y); // Y= (MdagM)^-1 Mdag Vdag phi = M^-1 V phi + NumOp.ExportFourDimPseudoFermion(Y,Y4); + + RealD action = norm2(Y4); + + return action; + }; + + ////////////////////////////////////////////////////// + // dS/du = 2 Re phi^dag (V^dag M^-dag)_11 (M^-1 d V)_11 phi + // - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + FermionField phi(NumOp.FermionGrid()); + FermionField Vphi(NumOp.FermionGrid()); + FermionField MinvVphi(NumOp.FermionGrid()); + FermionField tmp4(NumOp.GaugeGrid()); + FermionField MdagInvMinvVphi(NumOp.FermionGrid()); + + GaugeField force(NumOp.GaugeGrid()); + + //Y=V phi + //X = (Mdag V phi + //Y = (Mdag M)^-1 Mdag V phi = M^-1 V Phi + NumOp.ImportFourDimPseudoFermion(phi4,phi); + NumOp.M(phi,Vphi); // V phi + SchurRedBlackDiagMooeeSolve PrecSolve(DerivativeSolver); + PrecSolve(DenOp,Vphi,MinvVphi);// M^-1 V phi + std::cout << GridLogMessage << "4d deriv solve "<< norm2(MinvVphi)<<"\n"; + + // Projects onto the physical space and back + NumOp.ExportFourDimPseudoFermion(MinvVphi,tmp4); + NumOp.ImportFourDimPseudoFermion(tmp4,Y); + + SchurRedBlackDiagMooeeDagSolve PrecDagSolve(DerivativeDagSolver); + // X = proj M^-dag V phi + // Need an adjoint solve + PrecDagSolve(DenOp,Y,MdagInvMinvVphi); + std::cout << GridLogMessage << "4d deriv solve dag "<< norm2(MdagInvMinvVphi)<<"\n"; + + // phi^dag (Vdag Mdag^-1) (M^-1 dV) phi + NumOp.MDeriv(force ,MdagInvMinvVphi , phi, DaggerNo ); dSdU=force; + + // phi^dag (dVdag Mdag^-1) (M^-1 V) phi + NumOp.MDeriv(force , phi, MdagInvMinvVphi ,DaggerYes ); dSdU=dSdU+force; + + // - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi + DenOp.MDeriv(force,MdagInvMinvVphi,MinvVphi,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,MinvVphi,MdagInvMinvVphi,DaggerYes); dSdU=dSdU-force; + + dSdU *= -1.0; + //dSdU = - Ta(dSdU); + + }; +}; + +NAMESPACE_END(Grid); + + From f8d7d2389326d17c57c2333ddff455bda964b701 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:19:53 +0200 Subject: [PATCH 13/24] 4D pseudofermion options --- HMC/Mobius2p1fEOFA_4dPseudoFermion.cc | 473 ++++++++++++++++++ HMC/Mobius2p1fRHMC_4dPseudoFermion.cc | 197 ++++++++ ...bius2p1fRHMC_4dPseudoFermionSchurSolver.cc | 198 ++++++++ 3 files changed, 868 insertions(+) create mode 100644 HMC/Mobius2p1fEOFA_4dPseudoFermion.cc create mode 100644 HMC/Mobius2p1fRHMC_4dPseudoFermion.cc create mode 100644 HMC/Mobius2p1fRHMC_4dPseudoFermionSchurSolver.cc diff --git a/HMC/Mobius2p1fEOFA_4dPseudoFermion.cc b/HMC/Mobius2p1fEOFA_4dPseudoFermion.cc new file mode 100644 index 00000000..5b8731ae --- /dev/null +++ b/HMC/Mobius2p1fEOFA_4dPseudoFermion.cc @@ -0,0 +1,473 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: + +Copyright (C) 2015-2016 + +Author: Peter Boyle +Author: Guido Cossu +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 */ +#include +#include +#include + +#ifdef GRID_DEFAULT_PRECISION_DOUBLE +#define MIXED_PRECISION +#endif + +NAMESPACE_BEGIN(Grid); + + /* + * Need a plan for gauge field update for mixed precision in HMC (2x speed up) + * -- Store the single prec action operator. + * -- Clone the gauge field from the operator function argument. + * -- Build the mixed precision operator dynamically from the passed operator and single prec clone. + */ + + template + class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction { + public: + typedef typename FermionOperatorD::FermionField FieldD; + typedef typename FermionOperatorF::FermionField FieldF; + + using OperatorFunction::operator(); + + RealD Tolerance; + RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed + Integer MaxInnerIterations; + Integer MaxOuterIterations; + GridBase* SinglePrecGrid4; //Grid for single-precision fields + GridBase* SinglePrecGrid5; //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 + + FermionOperatorF &FermOpF; + FermionOperatorD &FermOpD;; + SchurOperatorF &LinOpF; + SchurOperatorD &LinOpD; + + Integer TotalInnerIterations; //Number of inner CG iterations + Integer TotalOuterIterations; //Number of restarts + Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step + + MixedPrecisionConjugateGradientOperatorFunction(RealD tol, + Integer maxinnerit, + Integer maxouterit, + GridBase* _sp_grid4, + GridBase* _sp_grid5, + FermionOperatorF &_FermOpF, + FermionOperatorD &_FermOpD, + SchurOperatorF &_LinOpF, + SchurOperatorD &_LinOpD): + LinOpF(_LinOpF), + LinOpD(_LinOpD), + FermOpF(_FermOpF), + FermOpD(_FermOpD), + Tolerance(tol), + InnerTolerance(tol), + MaxInnerIterations(maxinnerit), + MaxOuterIterations(maxouterit), + SinglePrecGrid4(_sp_grid4), + SinglePrecGrid5(_sp_grid5), + OuterLoopNormMult(100.) + { + /* Debugging instances of objects; references are stored + std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " < &LinOpU, const FieldD &src, FieldD &psi) { + + std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<(&LinOpU); + + // std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <_Mat)<_Mat)==&(LinOpD._Mat)); + + //////////////////////////////////////////////////////////////////////////////////// + // Must snarf a single precision copy of the gauge field in Linop_d argument + //////////////////////////////////////////////////////////////////////////////////// + typedef typename FermionOperatorF::GaugeField GaugeFieldF; + typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF; + typedef typename FermionOperatorD::GaugeField GaugeFieldD; + typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD; + + GridBase * GridPtrF = SinglePrecGrid4; + GridBase * GridPtrD = FermOpD.Umu.Grid(); + GaugeFieldF U_f (GridPtrF); + GaugeLinkFieldF Umu_f(GridPtrF); + // std::cout << " Dim gauge field "<Nd()<Nd()<(FermOpD.Umu, mu); + precisionChange(Umu_f,Umu_d); + PokeIndex(FermOpF.Umu, Umu_f, mu); + } + pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu); + pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu); + + //////////////////////////////////////////////////////////////////////////////////// + // Could test to make sure that LinOpF and LinOpD agree to single prec? + //////////////////////////////////////////////////////////////////////////////////// + /* + GridBase *Fgrid = psi._grid; + FieldD tmp2(Fgrid); + FieldD tmp1(Fgrid); + LinOpU.Op(src,tmp1); + LinOpD.Op(src,tmp2); + std::cout << " Double gauge field "<< norm2(FermOpD.Umu)< MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD); + std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" < HMCWrapper; + // MD.name = std::string("Leap Frog"); + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Force Gradient"); + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 12; + MD.trajL = 1.0; + + HMCparameters HMCparams; + HMCparams.StartTrajectory = 211; + HMCparams.Trajectories = 1000; + HMCparams.NoMetropolisUntil= 0; + // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.MD = MD; + HMCWrapper TheHMC(HMCparams); + + // Grid from the command line arguments --grid and --mpi + TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition + + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_EOFA4D_lat"; + CPparams.rng_prefix = "ckpoint_EOFA4D_rng"; + CPparams.saveInterval = 1; + CPparams.format = "IEEE64BIG"; + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + const int Ls = 16; + Real beta = 2.13; + Real light_mass = 0.01; + Real strange_mass = 0.04; + Real pv_mass = 1.0; + RealD M5 = 1.8; + RealD b = 1.0; + RealD c = 0.0; + + std::vector hasenbusch({ 0.1, 0.3, 0.6 }); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); + auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); + + Coordinate latt = GridDefaultLatt(); + Coordinate mpi = GridDefaultMpi(); + Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd()); + Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd()); + auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi); + auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF); + auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF); + auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF); + + IwasakiGaugeActionR GaugeAction(beta); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + LatticeGaugeFieldF UF(GridPtrF); + + // These lines are unecessary if BC are all periodic + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + FermionActionF::ImplParams ParamsF(boundary); + + double ActionStoppingCondition = 1e-10; + double DerivativeStoppingCondition = 1e-8; + double MaxCGIterations = 30000; + + //////////////////////////////////// + // Collect actions + //////////////////////////////////// + ActionLevel Level1(1); + ActionLevel Level2(8); + + //////////////////////////////////// + // Strange action + //////////////////////////////////// + typedef SchurDiagMooeeOperator LinearOperatorF; + typedef SchurDiagMooeeOperator LinearOperatorD; + + typedef SchurDiagMooeeDagOperator LinearOperatorDagF; + typedef SchurDiagMooeeDagOperator LinearOperatorDagD; + + typedef SchurDiagMooeeOperator LinearOperatorEOFAF; + typedef SchurDiagMooeeOperator LinearOperatorEOFAD; + + typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG; + typedef MixedPrecisionConjugateGradientOperatorFunction MxDagPCG; + typedef MixedPrecisionConjugateGradientOperatorFunction MxPCG_EOFA; + + // DJM: setup for EOFA ratio (Mobius) + OneFlavourRationalParams OFRp; + OFRp.lo = 0.1; + OFRp.hi = 25.0; + OFRp.MaxIter = 10000; + OFRp.tolerance= 1.0e-9; + OFRp.degree = 14; + OFRp.precision= 50; + + + MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); + MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c); + MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); + MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c); + + ConjugateGradient ActionCG(ActionStoppingCondition,MaxCGIterations); + ConjugateGradient DerivativeCG(DerivativeStoppingCondition,MaxCGIterations); +#ifdef MIXED_PRECISION + const int MX_inner = 1000; + // Mixed precision EOFA + LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L); + LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R); + LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF); + LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF); + + MxPCG_EOFA ActionCGL(ActionStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + Strange_Op_LF,Strange_Op_L, + Strange_LinOp_LF,Strange_LinOp_L); + + MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + Strange_Op_LF,Strange_Op_L, + Strange_LinOp_LF,Strange_LinOp_L); + + MxPCG_EOFA ActionCGR(ActionStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + Strange_Op_RF,Strange_Op_R, + Strange_LinOp_RF,Strange_LinOp_R); + + MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition, + MX_inner, + MaxCGIterations, + GridPtrF, + FrbGridF, + Strange_Op_RF,Strange_Op_R, + Strange_LinOp_RF,Strange_LinOp_R); + + ExactOneFlavourRatioPseudoFermionAction + EOFA(Strange_Op_L, Strange_Op_R, + ActionCG, + ActionCGL, ActionCGR, + DerivativeCGL, DerivativeCGR, + OFRp, true); +#else + ExactOneFlavourRatioPseudoFermionAction + EOFA(Strange_Op_L, Strange_Op_R, + ActionCG, + ActionCG, ActionCG, + DerivativeCG, DerivativeCG, + OFRp, true); +#endif + Level1.push_back(&EOFA); + + //////////////////////////////////// + // up down action + //////////////////////////////////// + std::vector light_den; + std::vector light_num; + + int n_hasenbusch = hasenbusch.size(); + light_den.push_back(light_mass); + for(int h=0;h Numerators; + std::vector Denominators; + std::vector *> Quotients; + std::vector ActionMPCG; + std::vector MPCG; + std::vector MPCGdag; + std::vector DenominatorsF; + std::vector LinOpD; + std::vector LinOpF; + std::vector LinOpDagD; + std::vector LinOpDagF; + + for(int h=0;h(*Numerators[h],*Denominators[h],*MPCG[h],*MPCGdag[h],*ActionMPCG[h],ActionCG)); +#else + //////////////////////////////////////////////////////////////////////////// + // Standard CG for 2f force + //////////////////////////////////////////////////////////////////////////// + Quotients.push_back (new TwoFlavourRatioEO4DPseudoFermionAction(*Numerators[h],*Denominators[h],DerivativeCG,ActionCG)); +#endif + + } + + for(int h=0;h +Author: Guido Cossu + +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 + +int main(int argc, char **argv) { + using namespace Grid; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef WilsonImplR FermionImplPolicy; + typedef MobiusFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + typedef Grid::XmlReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + IntegratorParameters MD; + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Leap Frog"); + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Force Gradient"); + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 10; + MD.trajL = 1.0; + + HMCparameters HMCparams; + HMCparams.StartTrajectory = 137; + HMCparams.Trajectories = 200; + HMCparams.NoMetropolisUntil= 0; + // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.MD = MD; + HMCWrapper TheHMC(HMCparams); + + // Grid from the command line arguments --grid and --mpi + TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition + + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_EODWF_lat"; + CPparams.rng_prefix = "ckpoint_EODWF_rng"; + CPparams.saveInterval = 1; + CPparams.format = "IEEE64BIG"; + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + const int Ls = 16; + Real beta = 2.13; + Real light_mass = 0.01; + Real strange_mass = 0.04; + Real pv_mass = 1.0; + RealD M5 = 1.8; + RealD b = 1.0; + RealD c = 0.0; + + // FIXME: + // Same in MC and MD + // Need to mix precision too + OneFlavourRationalParams OFRp; + OFRp.lo = 4.0e-3; + OFRp.hi = 30.0; + OFRp.MaxIter = 10000; + OFRp.tolerance= 1.0e-10; + OFRp.degree = 16; + OFRp.precision= 50; + + std::vector hasenbusch({ 0.1 }); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); + auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); + + IwasakiGaugeActionR GaugeAction(beta); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + // These lines are unecessary if BC are all periodic + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + + double StoppingCondition = 1e-10; + double MaxCGIterations = 30000; + ConjugateGradient CG(StoppingCondition,MaxCGIterations); + + //////////////////////////////////// + // Collect actions + //////////////////////////////////// + ActionLevel Level1(1); + ActionLevel Level2(8); + + //////////////////////////////////// + // Strange action + //////////////////////////////////// + + // FermionAction StrangeOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_mass,M5,b,c, Params); + // DomainWallEOFAFermionR Strange_Op_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5); + // DomainWallEOFAFermionR Strange_Op_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5); + // ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L,Strange_Op_R,CG,ofp, false); + + FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params); + FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params); + + OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermion(StrangePauliVillarsOp,StrangeOp,OFRp); + Level1.push_back(&StrangePseudoFermion); + + //////////////////////////////////// + // up down action + //////////////////////////////////// + std::vector light_den; + std::vector light_num; + + int n_hasenbusch = hasenbusch.size(); + light_den.push_back(light_mass); + for(int h=0;h Numerators; + std::vector Denominators; + std::vector *> Quotients; + + for(int h=0;h(*Numerators[h],*Denominators[h],CG,CG)); + } + + for(int h=0;h + +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 +#include +#include + +int main(int argc, char **argv) { + using namespace Grid; + + Grid_init(&argc, &argv); + int threads = GridThread::GetThreads(); + // here make a routine to print all the relevant information on the run + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + + // Typedefs to simplify notation + typedef WilsonImplR FermionImplPolicy; + typedef MobiusFermionR FermionAction; + typedef typename FermionAction::FermionField FermionField; + + typedef Grid::XmlReader Serialiser; + + //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + IntegratorParameters MD; + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Leap Frog"); + // typedef GenericHMCRunner HMCWrapper; + // MD.name = std::string("Force Gradient"); + typedef GenericHMCRunner HMCWrapper; + MD.name = std::string("MinimumNorm2"); + MD.MDsteps = 12; + MD.trajL = 1.0; + + HMCparameters HMCparams; + HMCparams.StartTrajectory = 211; + HMCparams.Trajectories = 200; + HMCparams.NoMetropolisUntil= 0; + // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; + // HMCparams.StartingType =std::string("ColdStart"); + HMCparams.StartingType =std::string("CheckpointStart"); + HMCparams.MD = MD; + HMCWrapper TheHMC(HMCparams); + + // Grid from the command line arguments --grid and --mpi + TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition + + CheckpointerParameters CPparams; + CPparams.config_prefix = "ckpoint_4dDWF_lat"; + CPparams.rng_prefix = "ckpoint_4dDWF_rng"; + CPparams.saveInterval = 1; + CPparams.format = "IEEE64BIG"; + TheHMC.Resources.LoadNerscCheckpointer(CPparams); + + RNGModuleParameters RNGpar; + RNGpar.serial_seeds = "1 2 3 4 5"; + RNGpar.parallel_seeds = "6 7 8 9 10"; + TheHMC.Resources.SetRNGSeeds(RNGpar); + + // Construct observables + // here there is too much indirection + typedef PlaquetteMod PlaqObs; + TheHMC.Resources.AddObservable(); + ////////////////////////////////////////////// + + const int Ls = 16; + Real beta = 2.13; + Real light_mass = 0.01; + Real strange_mass = 0.04; + Real pv_mass = 1.0; + RealD M5 = 1.8; + RealD b = 1.0; + RealD c = 0.0; + + // FIXME: + // Same in MC and MD + // Need to mix precision too + OneFlavourRationalParams OFRp; + OFRp.lo = 4.0e-3; + OFRp.hi = 30.0; + OFRp.MaxIter = 10000; + OFRp.tolerance= 1.0e-10; + OFRp.degree = 16; + OFRp.precision= 50; + + std::vector hasenbusch({ 0.1 }); + + auto GridPtr = TheHMC.Resources.GetCartesian(); + auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); + auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr); + auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr); + + IwasakiGaugeActionR GaugeAction(beta); + + // temporarily need a gauge field + LatticeGaugeField U(GridPtr); + + // These lines are unecessary if BC are all periodic + std::vector boundary = {1,1,1,-1}; + FermionAction::ImplParams Params(boundary); + + double StoppingCondition = 1e-10; + double MaxCGIterations = 30000; + ConjugateGradient CG(StoppingCondition,MaxCGIterations); + + //////////////////////////////////// + // Collect actions + //////////////////////////////////// + ActionLevel Level1(1); + ActionLevel Level2(8); + + //////////////////////////////////// + // Strange action + //////////////////////////////////// + + // FermionAction StrangeOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_mass,M5,b,c, Params); + // DomainWallEOFAFermionR Strange_Op_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5); + // DomainWallEOFAFermionR Strange_Op_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5); + // ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L,Strange_Op_R,CG,ofp, false); + + FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params); + FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params); + + OneFlavourEvenOddRatioRationalPseudoFermionAction StrangePseudoFermion(StrangePauliVillarsOp,StrangeOp,OFRp); + Level1.push_back(&StrangePseudoFermion); + + //////////////////////////////////// + // up down action + //////////////////////////////////// + std::vector light_den; + std::vector light_num; + + int n_hasenbusch = hasenbusch.size(); + light_den.push_back(light_mass); + for(int h=0;h Numerators; + std::vector Denominators; + std::vector *> Quotients; + + for(int h=0;h(*Numerators[h],*Denominators[h],CG,CG)); + } + + for(int h=0;h Date: Thu, 6 May 2021 23:21:15 +0200 Subject: [PATCH 14/24] Schur solver for Mdag --- Grid/algorithms/iterative/SchurRedBlack.h | 154 +++++++++++++++++++++- 1 file changed, 147 insertions(+), 7 deletions(-) diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h index d0b133a3..4706056a 100644 --- a/Grid/algorithms/iterative/SchurRedBlack.h +++ b/Grid/algorithms/iterative/SchurRedBlack.h @@ -40,7 +40,7 @@ Author: Peter Boyle * (-MoeMee^{-1} 1 ) * L^{dag} = ( 1 Mee^{-dag} Moe^{dag} ) * ( 0 1 ) - * L^{-d} = ( 1 -Mee^{-dag} Moe^{dag} ) + * L^{-dag}= ( 1 -Mee^{-dag} Moe^{dag} ) * ( 0 1 ) * * U^-1 = (1 -Mee^{-1} Meo) @@ -82,7 +82,8 @@ Author: Peter Boyle * c) M_oo^-dag Doo^{dag} Doo Moo^-1 phi_0 = M_oo^-dag (D_oo)^dag L^{-1} eta_o * eta_o' = M_oo^-dag (D_oo)^dag (eta_o - Moe Mee^{-1} eta_e) * psi_o = M_oo^-1 phi_o - * TODO: Deflation + * + * */ namespace Grid { @@ -97,6 +98,7 @@ namespace Grid { protected: typedef CheckerBoardedSparseMatrixBase Matrix; OperatorFunction & _HermitianRBSolver; + int CBfactorise; bool subGuess; bool useSolnAsInitGuess; // if true user-supplied solution vector is used as initial guess for solver @@ -189,13 +191,20 @@ namespace Grid { ///////////////////////////////////////////////// // Check unprec residual if possible ///////////////////////////////////////////////// - if ( ! subGuess ) { - _Matrix.M(out[b],resid); + if ( ! subGuess ) { + + if ( this->adjoint() ) _Matrix.Mdag(out[b],resid); + else _Matrix.M(out[b],resid); + resid = resid-in[b]; RealD ns = norm2(in[b]); RealD nr = norm2(resid); - std::cout<adjoint() << std::endl; + if ( this->adjoint() ) + std::cout<adjoint() << std::endl; + + if ( this->adjoint() ) _Matrix.Mdag(out,resid); + else _Matrix.M(out,resid); + resid = resid-in; RealD ns = norm2(in); RealD nr = norm2(resid); - std::cout<adjoint() ) + std::cout<_HermitianRBSolver(_OpEO, src_o, sol_o); } }; + + /* + * Red black Schur decomposition + * + * M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo) + * (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 ) + * = L D U + * + * L^-1 = (1 0 ) + * (-MoeMee^{-1} 1 ) + * L^{dag} = ( 1 Mee^{-dag} Moe^{dag} ) + * ( 0 1 ) + * + * U^-1 = (1 -Mee^{-1} Meo) + * (0 1 ) + * U^{dag} = ( 1 0) + * (Meo^dag Mee^{-dag} 1) + * U^{-dag} = ( 1 0) + * (-Meo^dag Mee^{-dag} 1) + * + * + *********************** + * M^dag psi = eta + *********************** + * + * Really for Mobius: (Wilson - easier to just use gamma 5 hermiticity) + * + * Mdag psi = Udag Ddag Ldag psi = eta + * + * U^{-dag} = ( 1 0) + * (-Meo^dag Mee^{-dag} 1) + * + * + * i) D^dag phi = (U^{-dag} eta) + * eta'_e = eta_e + * eta'_o = (eta_o - Meo^dag Mee^{-dag} eta_e) + * + * phi_o = D_oo^-dag eta'_o = D_oo^-dag (eta_o - Meo^dag Mee^{-dag} eta_e) + * + * phi_e = D_ee^-dag eta'_e = D_ee^-dag eta_e + * + * Solve: + * + * D_oo D_oo^dag phi_o = D_oo (eta_o - Meo^dag Mee^{-dag} eta_e) + * + * ii) + * phi = L^dag psi => psi = L^-dag phi. + * + * L^{-dag} = ( 1 -Mee^{-dag} Moe^{dag} ) + * ( 0 1 ) + * + * => sol_e = M_ee^-dag * ( src_e - Moe^dag phi_o )... + * => sol_o = phi_o + */ + /////////////////////////////////////////////////////////////////////////////////////////////////////// + // Site diagonal has Mooee on it, but solve the Adjoint system + /////////////////////////////////////////////////////////////////////////////////////////////////////// + template class SchurRedBlackDiagMooeeDagSolve : public SchurRedBlackBase { + public: + typedef CheckerBoardedSparseMatrixBase Matrix; + + virtual bool adjoint(void) { return true; } + SchurRedBlackDiagMooeeDagSolve(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 = (source_o - Moe^dag MeeInvDag source_e) + ///////////////////////////////////////////////////// + _Matrix.MooeeInvDag(src_e,tmp); assert( tmp.Checkerboard() ==Even); + _Matrix.MeooeDag (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd); + tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd); + + // get the right Mpc + SchurDiagMooeeOperator _HermOpEO(_Matrix); + _HermOpEO.Mpc(tmp,src_o); assert(src_o.Checkerboard() ==Odd); + } + virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) + { + SchurDiagMooeeDagOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + }; + virtual void RedBlackSolve (Matrix & _Matrix,const std::vector &src_o, std::vector &sol_o) + { + SchurDiagMooeeDagOperator _HermOpEO(_Matrix); + this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); + } + 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_e(grid); + Field tmp(grid); + /////////////////////////////////////////////////// + // sol_e = M_ee^-dag * ( src_e - Moe^dag phi_o )... + // sol_o = phi_o + /////////////////////////////////////////////////// + _Matrix.MeooeDag(sol_o,tmp); assert(tmp.Checkerboard()==Even); + tmp = src_e-tmp; assert(tmp.Checkerboard()==Even); + _Matrix.MooeeInvDag(tmp,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 ); + } + }; + } #endif From 3b433fe6fba7ada65160b326021b043b56b0ad7b Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:22:09 +0200 Subject: [PATCH 15/24] Better force logging --- Grid/qcd/hmc/integrators/Integrator.h | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h index 6b0e3caf..27b66fae 100644 --- a/Grid/qcd/hmc/integrators/Integrator.h +++ b/Grid/qcd/hmc/integrators/Integrator.h @@ -132,11 +132,20 @@ protected: as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl; + auto name = as[level].actions.at(a)->action_name(); if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); force = FieldImplementation::projectForce(force); // Ta for gauge fields double end_force = usecond(); Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites()); - std::cout << GridLogIntegrator << "["<gSites()); + std::cout << GridLogIntegrator + << "["< Date: Thu, 6 May 2021 23:23:30 +0200 Subject: [PATCH 16/24] Code prettier --- Grid/qcd/utils/CovariantCshift.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h index 6c70706f..29f796a4 100644 --- a/Grid/qcd/utils/CovariantCshift.h +++ b/Grid/qcd/utils/CovariantCshift.h @@ -182,7 +182,7 @@ namespace ConjugateBC { GridBase *grid = Link.Grid(); int Lmu = grid->GlobalDimensions()[mu] - 1; - Lattice> coor(grid); + Lattice > coor(grid); LatticeCoordinate(coor, mu); Lattice tmp(grid); From a0534e03f913149166b103e7e583dfc0452f1e96 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:23:57 +0200 Subject: [PATCH 17/24] Augmented test --- tests/debug/Test_cayley_even_odd.cc | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/debug/Test_cayley_even_odd.cc b/tests/debug/Test_cayley_even_odd.cc index 5e800b26..5688894e 100644 --- a/tests/debug/Test_cayley_even_odd.cc +++ b/tests/debug/Test_cayley_even_odd.cc @@ -208,6 +208,27 @@ void TestWhat(What & Ddwf, std::cout< * = < chi | Deo^dag| phi> "< Date: Thu, 6 May 2021 23:24:19 +0200 Subject: [PATCH 18/24] Fix compile --- tests/debug/Test_heatbath_dwf_eofa.cc | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/tests/debug/Test_heatbath_dwf_eofa.cc b/tests/debug/Test_heatbath_dwf_eofa.cc index 9d453a96..d3b3cbe9 100644 --- a/tests/debug/Test_heatbath_dwf_eofa.cc +++ b/tests/debug/Test_heatbath_dwf_eofa.cc @@ -66,10 +66,9 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); - GridParallelRNG RNG5(FGrid); - RNG5.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4(UGrid); - RNG4.SeedFixedIntegers(seeds4); + GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); // Random gauge field LatticeGaugeField Umu(UGrid); @@ -84,7 +83,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -94,7 +93,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } From 7533f66b5457bea6c94d8d535efe26ec400a1fa7 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:24:39 +0200 Subject: [PATCH 19/24] Fix compile --- tests/debug/Test_heatbath_dwf_eofa_gparity.cc | 11 +++++------ tests/debug/Test_heatbath_mobius_eofa.cc | 5 +++-- tests/debug/Test_heatbath_mobius_eofa_gparity.cc | 5 +++-- 3 files changed, 11 insertions(+), 10 deletions(-) diff --git a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc index 22cc1e90..123553a6 100644 --- a/tests/debug/Test_heatbath_dwf_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_dwf_eofa_gparity.cc @@ -70,10 +70,9 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); - GridParallelRNG RNG5(FGrid); - RNG5.SeedFixedIntegers(seeds5); - GridParallelRNG RNG4(UGrid); - RNG4.SeedFixedIntegers(seeds4); + GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); // Random gauge field LatticeGaugeField Umu(UGrid); @@ -90,7 +89,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -100,7 +99,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa.cc b/tests/debug/Test_heatbath_mobius_eofa.cc index 4cf4bf53..b4b7b983 100644 --- a/tests/debug/Test_heatbath_mobius_eofa.cc +++ b/tests/debug/Test_heatbath_mobius_eofa.cc @@ -68,6 +68,7 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds4); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); @@ -86,7 +87,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -96,7 +97,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } diff --git a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc index 2fcb4b9f..ff5ec81c 100644 --- a/tests/debug/Test_heatbath_mobius_eofa_gparity.cc +++ b/tests/debug/Test_heatbath_mobius_eofa_gparity.cc @@ -72,6 +72,7 @@ int main(int argc, char** argv) // Set up RNGs std::vector seeds4({1, 2, 3, 4}); std::vector seeds5({5, 6, 7, 8}); + GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); @@ -91,7 +92,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, false); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } @@ -101,7 +102,7 @@ int main(int argc, char** argv) ConjugateGradient CG(1.0e-12, 5000); ExactOneFlavourRatioPseudoFermionAction Meofa(Lop, Rop, CG, Params, true); - Meofa.refresh(Umu, RNG5); + Meofa.refresh(Umu, SRNG, RNG5); printf(" = %1.15e\n", Meofa.S(Umu)); } From bab88bc4f7bd59bbeb4e088e3bd933681e82c2e1 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:25:12 +0200 Subject: [PATCH 20/24] Fix compile --- tests/forces/Test_gp_rect_force.cc | 11 ++++++----- tests/forces/Test_momentum_filter.cc | 5 ++++- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/tests/forces/Test_gp_rect_force.cc b/tests/forces/Test_gp_rect_force.cc index 98ebb2fa..627d264e 100644 --- a/tests/forces/Test_gp_rect_force.cc +++ b/tests/forces/Test_gp_rect_force.cc @@ -29,7 +29,6 @@ Author: paboyle using namespace std; using namespace Grid; - ; @@ -53,12 +52,14 @@ int main (int argc, char ** argv) pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); LatticeGaugeField U(&Grid); - SU::HotConfiguration(pRNG,U); - + double beta = 1.0; double c1 = 0.331; + std::vector twists(Nd,0); + twists[1] = 0; + ConjugateGimplD::setDirections(twists); ConjugatePlaqPlusRectangleActionR Action(beta,c1); //ConjugateWilsonGaugeActionR Action(beta); //WilsonGaugeActionR Action(beta); @@ -67,13 +68,13 @@ int main (int argc, char ** argv) // get the deriv of phidag MdagM phi with respect to "U" LatticeGaugeField UdSdU(&Grid); - + Action.deriv(U,UdSdU); //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.01; + RealD dt = 0.001; LatticeColourMatrix mommu(&Grid); LatticeColourMatrix forcemu(&Grid); diff --git a/tests/forces/Test_momentum_filter.cc b/tests/forces/Test_momentum_filter.cc index 856ea0f2..c7888948 100644 --- a/tests/forces/Test_momentum_filter.cc +++ b/tests/forces/Test_momentum_filter.cc @@ -59,9 +59,12 @@ int main (int argc, char ** argv) std::cout< seeds({1,2,3,4}); + std::vector serial_seeds({5,6,7,8}); + GridSerialRNG sRNG; GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds); + sRNG.SeedFixedIntegers(serial_seeds); typedef PeriodicGimplR Gimpl; typedef WilsonGaugeAction GaugeAction; @@ -115,7 +118,7 @@ int main (int argc, char ** argv) integrator.setMomentumFilter(filter); - integrator.refresh(U, pRNG); //doesn't actually change the gauge field + integrator.refresh(U, sRNG, pRNG); //doesn't actually change the gauge field //Check the momentum is zero on the boundary const auto &P = integrator.getMomentum(); From 374fb325f3bb514bbfe8f1cae44fdf49aab41688 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:26:42 +0200 Subject: [PATCH 21/24] Remove LIME dependency --- tests/hmc/Test_hmc_GparityIwasakiGauge.cc | 8 +++++--- tests/hmc/Test_hmc_IwasakiGauge.cc | 6 ++++-- tests/hmc/Test_hmc_ScalarActionNxN.cc | 4 ++-- tests/hmc/Test_hmc_WG_Production.cc | 4 ++-- 4 files changed, 13 insertions(+), 9 deletions(-) diff --git a/tests/hmc/Test_hmc_GparityIwasakiGauge.cc b/tests/hmc/Test_hmc_GparityIwasakiGauge.cc index 7f74d5d8..c396eb1b 100644 --- a/tests/hmc/Test_hmc_GparityIwasakiGauge.cc +++ b/tests/hmc/Test_hmc_GparityIwasakiGauge.cc @@ -58,7 +58,7 @@ int main(int argc, char **argv) { CheckpointerParameters CPparams; CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.rng_prefix = "ckpoint_EODWF_rng"; - CPparams.saveInterval = 5; + CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); @@ -79,7 +79,7 @@ int main(int argc, char **argv) { // that have a complex construction // standard RealD beta = 2.6 ; - const int nu = 3; + const int nu = 1; std::vector twists(Nd,0); twists[nu] = 1; ConjugateGimplD::setDirections(twists); @@ -95,8 +95,10 @@ int main(int argc, char **argv) { ///////////////////////////////////////////////////////////// // HMC parameters are serialisable - TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.MDsteps = 80; TheHMC.Parameters.MD.trajL = 1.0; + TheHMC.Parameters.Trajectories = 100; + TheHMC.Parameters.NoMetropolisUntil= 10; TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file diff --git a/tests/hmc/Test_hmc_IwasakiGauge.cc b/tests/hmc/Test_hmc_IwasakiGauge.cc index 66345bc0..a5484280 100644 --- a/tests/hmc/Test_hmc_IwasakiGauge.cc +++ b/tests/hmc/Test_hmc_IwasakiGauge.cc @@ -51,7 +51,7 @@ int main(int argc, char **argv) { CheckpointerParameters CPparams; CPparams.config_prefix = "ckpoint_lat"; CPparams.rng_prefix = "ckpoint_rng"; - CPparams.saveInterval = 20; + CPparams.saveInterval = 1; CPparams.format = "IEEE64BIG"; TheHMC.Resources.LoadNerscCheckpointer(CPparams); @@ -81,8 +81,10 @@ int main(int argc, char **argv) { ///////////////////////////////////////////////////////////// // HMC parameters are serialisable - TheHMC.Parameters.MD.MDsteps = 20; + TheHMC.Parameters.MD.MDsteps = 80; TheHMC.Parameters.MD.trajL = 1.0; + TheHMC.Parameters.Trajectories = 100; + TheHMC.Parameters.NoMetropolisUntil= 10; TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file TheHMC.Run(); // no smearing diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index 0c398306..d2bf2955 100644 --- a/tests/hmc/Test_hmc_ScalarActionNxN.cc +++ b/tests/hmc/Test_hmc_ScalarActionNxN.cc @@ -132,8 +132,8 @@ int main(int argc, char **argv) { // Checkpointer definition CheckpointerParameters CPparams(Reader); - //TheHMC.Resources.LoadBinaryCheckpointer(CPparams); - TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar); + TheHMC.Resources.LoadBinaryCheckpointer(CPparams); + //TheHMC.Resources.LoadNerscCheckpointer(CPparams, SPar); RNGModuleParameters RNGpar(Reader); TheHMC.Resources.SetRNGSeeds(RNGpar); diff --git a/tests/hmc/Test_hmc_WG_Production.cc b/tests/hmc/Test_hmc_WG_Production.cc index b16c073b..aa731034 100644 --- a/tests/hmc/Test_hmc_WG_Production.cc +++ b/tests/hmc/Test_hmc_WG_Production.cc @@ -74,10 +74,10 @@ int main(int argc, char **argv) { // Checkpointer definition CheckpointerParameters CPparams(Reader); - //TheHMC.Resources.LoadNerscCheckpointer(CPparams); + TheHMC.Resources.LoadNerscCheckpointer(CPparams); // Store metadata in the Scidac checkpointer - TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar); + //TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar); RNGModuleParameters RNGpar(Reader); TheHMC.Resources.SetRNGSeeds(RNGpar); From ffcab648907b0dccfbe8b603028a3de52b74c0d1 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:27:57 +0200 Subject: [PATCH 22/24] DDHMC filter --- Grid/qcd/action/momentum/DDHMCfilter.h | 79 ++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 Grid/qcd/action/momentum/DDHMCfilter.h diff --git a/Grid/qcd/action/momentum/DDHMCfilter.h b/Grid/qcd/action/momentum/DDHMCfilter.h new file mode 100644 index 00000000..27e0ac66 --- /dev/null +++ b/Grid/qcd/action/momentum/DDHMCfilter.h @@ -0,0 +1,79 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/hmc/DDHMC.h + +Copyright (C) 2021 + +Author: Peter Boyle +Author: Christopher Kelly + +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 */ + +//////////////////////////////////////////////////// +// DDHMC filter with sub-block size B[mu] +//////////////////////////////////////////////////// + +template +struct DDHMCFilter: public MomentumFilterBase +{ + typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename MomentaField::scalar_type scalar_type; //scalar complex type + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef iScalar > > ScalarType; //complex phase for each site + typedef Lattice LatticeLorentzScalarType; + typedef Lattice LatticeScalarType; + + Coordinate Block; + + // Could also zero links in plane like Luscher advocates. + + DDHMCFilter(const Coordinate &_Block): Block(_Block){} + + void applyFilter(MomentaField &P) const override + { + GridBase *grid = P.Grid(); + LatticeScalarType mask_mu(grid); + LatticeLorentzScalarType mask(grid); + + //////////////////////////////////////////////////// + // Zero strictly links crossing between domains + // Luscher also zeroes links in plane of domain boundaries + // Keeping interior only. This prevents force from plaquettes + // crossing domains and keeps whole MD trajectory local. + // Step further than I want to go. + //////////////////////////////////////////////////// + ScalarType zz = scalar_type(0.0); + ScalarType one= scalar_type(1.0); + + LatticeInteger coor(grid); + for(int mu=0;mu(mask, mask_mu, mu); + + } + } +}; + From b32fd473f819516f9193eac130548f86e58bd1d2 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:28:28 +0200 Subject: [PATCH 23/24] Start at the Domain decomposed supprt --- .../pseudofermion/DomainDecomposedBoundary.h | 346 ++++++++++++++++++ 1 file changed, 346 insertions(+) create mode 100644 Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h diff --git a/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h b/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h new file mode 100644 index 00000000..36a849e0 --- /dev/null +++ b/Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h @@ -0,0 +1,346 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h + + Copyright (C) 2015 + +Author: Peter Boyle +Author: Peter Boyle +Author: paboyle + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#pragma once + +NAMESPACE_BEGIN(Grid); + +/////////////////////////////////////// +// Two flavour ratio +/////////////////////////////////////// +template +class DomainDecomposedBoundary { +public: + INHERIT_IMPL_TYPES(Impl); + + typedef typename GaugeField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename GaugeField::scalar_type scalar_type; //scalar complex type + + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef iScalar > > ScalarType; //complex phase for each site + typedef Lattice LatticeLorentzScalarType; + typedef Lattice LatticeScalarType; + + Coordinate Block; + DDHMCFilter Filter; + const int Omega=0; + const int OmegaBar=1; + + void ProjectBoundaryBothDomains (FermionField &f,int sgn) + { + assert((sgn==1)||(sgn==-1)); + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + GridBase *grid = f.Grid(); + LatticeInteger coor(grid); + LatticeInteger face(grid); + LatticeInteger nface(grid); nface=Zero(); + + ComplexField zz(grid); zz=Zero(); + + FermionField projected(grid); projected=Zero(); + FermionField sp_proj (grid); + + int dims = grid->Nd(); + int isDWF= (dims==Nd+1); + assert((dims==Nd)||(dims==Nd+1)); + + for(int mu=0;mu1,f,projected); + } + void ProjectDomain(FermionField &f,int cb) + { + GridBase *grid = f.Grid(); + ComplexField zz(grid); zz=Zero(); + LatticeInteger coor(grid); + LatticeInteger domaincb(grid); domaincb=Zero(); + for(int d=0;dNd();d++){ + LatticeCoordinate(coor,mu); + domaincb = domaincb + div(coor,Block[d]); + } + f = where(mod(domaincb,2)==cb,f,zz); + }; + + void ProjectOmegaBar (FermionField &f) {ProjectDomain(f,OmegaBar);} + void ProjectOmega (FermionField &f) {ProjectDomain(f,Omega);} + + // See my notes(!). + // Notation: Following Luscher, we introduce projectors $\hPdb$ with both spinor and space structure + // projecting all spinor elements in $\Omega$ connected by $\Ddb$ to $\bar{\Omega}$, + void ProjectBoundaryBar(FermionField &f) + { + ProjectBoundaryBothDomains(f); + ProjectOmega(f); + } + // and $\hPd$ projecting all spinor elements in $\bar{\Omega}$ connected by $\Dd$ to $\Omega$. + void ProjectBoundary (FermionField &f) + { + ProjectBoundaryBothDomains(f); + ProjectOmegaBar(f); + }; + + void dBoundary (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmegaBar(tmp); + Op.M(tmp,out); + ProjectOmega(out); + }; + void dBoundaryBar (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmega(tmp); + Op.M(tmp,out); + ProjectOmegaBar(out); + }; + void dOmega (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmega(tmp); + Op.M(tmp,out); + ProjectOmega(out); + }; + void dOmegaBar (FermionOperator &Op,FermionField &in,FermionField &out) + { + FermionField tmp(in); + ProjectOmegaBar(tmp); + Op.M(tmp,out); + ProjectOmegaBar(out); + }; + + void SolveOmega (FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void SolveOmegaBar(FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void SolveOmegaAndOmegaBar(FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + void dInverse (FermionOperator &Op,FermionField &in,FermionField &out){ assert(0); }; + + // R = Pdbar - Pdbar DomegaInv Dd DomegabarInv Ddbar + void R(FermionOperator &Op,FermionOperator &OpDirichlet,FermionField &in,FermionField &out) + { + FermionField tmp1(Op.FermionGrid()); + FermionField tmp2(Op.FermionGrid()); + dBoundaryBar(Op,in,tmp1); + SolveOmegaBar(OpDirichlet,tmp1,tmp2); // 1/2 cost + dBoundary(Op,tmp2,tmp1); + SolveOmega(OpDirichlet,tmp1,tmp2); // 1/2 cost + out = in - tmp2 ; + ProjectBoundaryBar(out); + }; + + // R = Pdbar - Pdbar Dinv Ddbar + void Rinverse(FermionField &in,FermionField &out) + { + FermionField tmp1(NumOp.FermionGrid()); + out = in; + ProjectBoundaryBar(out); + dInverse(out,tmp1); + ProjectBoundaryBar(tmp1); + out = out -tmp1; + }; + +} + +template +class DomainDecomposedBoundaryPseudoFermionAction : public Action { +public: + INHERIT_IMPL_TYPES(Impl); + +private: + FermionOperator & NumOp;// the basic operator + FermionOperator & DenOp;// the basic operator + FermionOperator & NumOpDirichlet;// the basic operator + FermionOperator & DenOpDirichlet;// the basic operator + + OperatorFunction &DerivativeSolver; + OperatorFunction &ActionSolver; + + FermionField Phi; // the pseudo fermion field for this trajectory + + +public: + DomainBoundaryPseudoFermionAction(FermionOperator &_NumOp, + FermionOperator &_DenOp, + FermionOperator &_NumOpDirichlet, + FermionOperator &_DenOpDirichlet, + OperatorFunction & DS, + OperatorFunction & AS, + Coordinate &_Block + ) : NumOp(_NumOp), + DenOp(_DenOp), + DerivativeSolver(DS), + ActionSolver(AS), + Phi(_NumOp.FermionGrid()), + Block(_Block) + // LinkFilter(Block) + {}; + + virtual std::string action_name(){return "DomainBoundaryPseudoFermionRatioAction";} + + virtual std::string LogParameters(){ + std::stringstream sstream; + sstream << GridLogMessage << "["< sig^2 = 0.5. + // + // So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707.... + // + RealD scale = std::sqrt(0.5); + + FermionField eta(NumOp.FermionGrid()); + FermionField tmp(NumOp.FermionGrid()); + + gaussian(pRNG,eta); + + ProjectBoundary(eta); + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + // Note: this hard codes normal equations type solvers; alternate implementation needed for + // non-herm style solvers. + MdagMLinearOperator ,FermionField> MdagMOp(NumOp); + + DenOp.Mdag(eta,Phi); // Mdag eta + tmp = Zero(); + ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta + NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta + + Phi=Phi*scale; + + }; + + ////////////////////////////////////////////////////// + // S = phi^dag V (Mdag M)^-1 Vdag phi + ////////////////////////////////////////////////////// + virtual RealD S(const GaugeField &U) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=Zero(); + ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + + RealD action = norm2(Y); + + return action; + }; + + ////////////////////////////////////////////////////// + // dS/du = phi^dag dV (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi + // + phi^dag V (Mdag M)^-1 dV^dag phi + ////////////////////////////////////////////////////// + virtual void deriv(const GaugeField &U,GaugeField & dSdU) { + + NumOp.ImportGauge(U); + DenOp.ImportGauge(U); + + MdagMLinearOperator ,FermionField> MdagMOp(DenOp); + + FermionField X(NumOp.FermionGrid()); + FermionField Y(NumOp.FermionGrid()); + + GaugeField force(NumOp.GaugeGrid()); + + + //Y=Vdag phi + //X = (Mdag M)^-1 V^dag phi + //Y = (Mdag)^-1 V^dag phi + NumOp.Mdag(Phi,Y); // Y= Vdag phi + X=Zero(); + DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi + DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi + + // phi^dag V (Mdag M)^-1 dV^dag phi + NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force; + + // phi^dag dV (Mdag M)^-1 V^dag phi + NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force; + + // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi + // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi + DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force; + DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force; + + dSdU *= -1.0; + //dSdU = - Ta(dSdU); + + }; +}; + +NAMESPACE_END(Grid); + + From 919ced1c3158af86c94f976e8fae5a4a0ed8be59 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 6 May 2021 23:29:03 +0200 Subject: [PATCH 24/24] Move the momentum filter to a better location for DDHMC --- Grid/qcd/action/momentum/MomentumFilter.h | 94 +++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 Grid/qcd/action/momentum/MomentumFilter.h diff --git a/Grid/qcd/action/momentum/MomentumFilter.h b/Grid/qcd/action/momentum/MomentumFilter.h new file mode 100644 index 00000000..2a15d80c --- /dev/null +++ b/Grid/qcd/action/momentum/MomentumFilter.h @@ -0,0 +1,94 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./lib/qcd/hmc/integrators/MomentumFilter.h + +Copyright (C) 2015 + +Author: Christopher Kelly +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 */ +//-------------------------------------------------------------------- +#ifndef MOMENTUM_FILTER +#define MOMENTUM_FILTER + +NAMESPACE_BEGIN(Grid); + +//These filter objects allow the user to manipulate the conjugate momentum as part of the update / refresh + +template +struct MomentumFilterBase{ + virtual void applyFilter(MomentaField &P) const; +}; + +//Do nothing +template +struct MomentumFilterNone: public MomentumFilterBase{ + void applyFilter(MomentaField &P) const override{} +}; + +//Multiply each site/direction by a Lorentz vector complex number field +//Can be used to implement a mask, zeroing out sites +template +struct MomentumFilterApplyPhase: public MomentumFilterBase{ + typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type + typedef typename MomentaField::scalar_type scalar_type; //scalar complex type + typedef iVector >, Nd > LorentzScalarType; //complex phase for each site/direction + typedef Lattice LatticeLorentzScalarType; + + LatticeLorentzScalarType phase; + + MomentumFilterApplyPhase(const LatticeLorentzScalarType _phase): phase(_phase){} + + //Default to uniform field of (1,0) + MomentumFilterApplyPhase(GridBase* _grid): phase(_grid){ + LorentzScalarType one; + for(int mu=0;mu