diff --git a/Grid/qcd/action/fermion/CloverHelpers.h b/Grid/qcd/action/fermion/CloverHelpers.h new file mode 100644 index 00000000..432bdf3c --- /dev/null +++ b/Grid/qcd/action/fermion/CloverHelpers.h @@ -0,0 +1,371 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h + + Copyright (C) 2017 - 2022 + + Author: paboyle + Author: Daniel Richtmann + Author: Mattia Bruno + + 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 + +#include +#include +#include + +//////////////////////////////////////////// +// Standard Clover +// (4+m0) + csw * clover_term +// Exp Clover +// (4+m0) * exp(csw/(4+m0) clover_term) +// = (4+m0) + csw * clover_term + ... +//////////////////////////////////////////// + +NAMESPACE_BEGIN(Grid); + + +////////////////////////////////// +// Generic Standard Clover +////////////////////////////////// + +template +class CloverHelpers: public WilsonCloverHelpers { +public: + + INHERIT_IMPL_TYPES(Impl); + INHERIT_CLOVER_TYPES(Impl); + + typedef WilsonCloverHelpers Helpers; + + static void Instantiate(CloverField& CloverTerm, CloverField& CloverTermInv, RealD csw_t, RealD diag_mass) { + GridBase *grid = CloverTerm.Grid(); + CloverTerm += diag_mass; + + int lvol = grid->lSites(); + int DimRep = Impl::Dimension; + { + autoView(CTv,CloverTerm,CpuRead); + autoView(CTIv,CloverTermInv,CpuWrite); + thread_for(site, lvol, { + Coordinate lcoor; + grid->LocalIndexToLocalCoor(site, lcoor); + Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); + typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero(); + peekLocalSite(Qx, CTv, lcoor); + //if (csw!=0){ + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++){ + auto zz = Qx()(j, k)(a, b); + EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex(zz); + } + // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; + + EigenInvCloverOp = EigenCloverOp.inverse(); + //std::cout << EigenInvCloverOp << std::endl; + for (int j = 0; j < Ns; j++) + for (int k = 0; k < Ns; k++) + for (int a = 0; a < DimRep; a++) + for (int b = 0; b < DimRep; b++) + Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); + // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; + // } + pokeLocalSite(Qxinv, CTIv, lcoor); + }); + } + } + + + static void CloverTermDerivative(GaugeField& clover_force, + const FermionField& X, + const FermionField& Y, + const std::vector& U, + RealD csw_t, RealD csw_r) { + GaugeLinkField force_mu(clover_force.Grid()), lambda(clover_force.Grid()); + PropagatorField Lambda(clover_force.Grid()); + + /////////////////////////////////////////////////////////// + // Clover term derivative + /////////////////////////////////////////////////////////// + Impl::outerProductImpl(Lambda, X, Y); + //std::cout << "Lambda:" << Lambda << std::endl; + + Gamma::Algebra sigma[] = { + Gamma::Algebra::SigmaXY, + Gamma::Algebra::SigmaXZ, + Gamma::Algebra::SigmaXT, + Gamma::Algebra::MinusSigmaXY, + Gamma::Algebra::SigmaYZ, + Gamma::Algebra::SigmaYT, + Gamma::Algebra::MinusSigmaXZ, + Gamma::Algebra::MinusSigmaYZ, + Gamma::Algebra::SigmaZT, + Gamma::Algebra::MinusSigmaXT, + Gamma::Algebra::MinusSigmaYT, + Gamma::Algebra::MinusSigmaZT}; + + /* + sigma_{\mu \nu}= + | 0 sigma[0] sigma[1] sigma[2] | + | sigma[3] 0 sigma[4] sigma[5] | + | sigma[6] sigma[7] 0 sigma[8] | + | sigma[9] sigma[10] sigma[11] 0 | + */ + + int count = 0; + clover_force = Zero(); + for (int mu = 0; mu < 4; mu++) + { + force_mu = Zero(); + for (int nu = 0; nu < 4; nu++) + { + if (mu == nu) + continue; + + RealD factor; + if (nu == 4 || mu == 4) + { + factor = 2.0 * csw_t; + } + else + { + factor = 2.0 * csw_r; + } + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu); // checked + count++; + } + + pokeLorentz(clover_force, U[mu] * force_mu, mu); + } + } +}; + + +////////////////////////////////// +// Generic Exp Clover +////////////////////////////////// + +template +class ExpCloverHelpers: public WilsonCloverHelpers { +public: + + INHERIT_IMPL_TYPES(Impl); + INHERIT_CLOVER_TYPES(Impl); + + template using iImplClover = iScalar, Ns>>; + typedef WilsonCloverHelpers Helpers; + + static void plusIdentity(const CloverField& in) { + int DimRep = Impl::Dimension; + + autoView(in_v, in, AcceleratorWrite); + + accelerator_for(ss, in.Grid()->oSites(), 1, { + for (int sa=0; saprec) { + NMAX++; + cond*=R/(double)(NMAX+1); + } + return NMAX; + } + + static int getNMAX(Lattice> &t, RealD R) {return getNMAX(1e-12,R);} + static int getNMAX(Lattice> &t, RealD R) {return getNMAX(1e-6,R);} + + static void Instantiate(CloverField& Clover, CloverField& CloverInv, RealD csw_t, RealD diag_mass) { + GridBase* grid = Clover.Grid(); + CloverField ExpClover(grid); + + int NMAX = getNMAX(Clover, 3.*csw_t/diag_mass); + + // csw/(diag_mass) * clover + Clover *= (1.0/diag_mass); + + ExpClover = Clover; + plusIdentity(ExpClover); // 1 + Clover + + // Taylor expansion, slow but generic + RealD pref = 1.0; + for (int i=2; i<=NMAX; i++) { + Clover = Clover * Clover; + pref /= (RealD)(i); + ExpClover += pref * Clover; + } + + // Convert the data layout of the clover terms + Clover = ExpClover * diag_mass; + CloverInv = adj(ExpClover * (1.0/diag_mass)); + } + + static void CloverTermDerivative(GaugeField& clover_force, const FermionField& X, const FermionField& Y, + const std::vector& U, RealD csw_t, RealD csw_r) { + assert(0); + } +}; + + +////////////////////////////////// +// Compact Standard Clover +////////////////////////////////// + + +template +class CompactCloverHelpers: public CompactWilsonCloverHelpers, + public WilsonCloverHelpers { +public: + + INHERIT_IMPL_TYPES(Impl); + INHERIT_CLOVER_TYPES(Impl); + INHERIT_COMPACT_CLOVER_TYPES(Impl); + + typedef WilsonCloverHelpers Helpers; + typedef CompactWilsonCloverHelpers CompactHelpers; + + static void MassTerm(CloverField& Clover, RealD diag_mass) { + Clover += diag_mass; + } + + static void Instantiate(CloverDiagonalField& Diagonal, + CloverTriangleField& Triangle, + CloverDiagonalField& DiagonalInv, + CloverTriangleField& TriangleInv, + RealD csw_t, RealD diag_mass) { + // Invert the clover term in the improved layout + CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv); + } + + // TODO: implement Cmunu for better performances with compact layout, but don't do it + // here, but rather in WilsonCloverHelpers.h -> CompactWilsonCloverHelpers + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) { + return Helpers::Cmunu(U, lambda, mu, nu); + } +}; + +////////////////////////////////// +// Compact Exp Clover +////////////////////////////////// + + +template +class CompactExpCloverHelpers: public CompactWilsonCloverHelpers { +public: + + INHERIT_IMPL_TYPES(Impl); + INHERIT_CLOVER_TYPES(Impl); + INHERIT_COMPACT_CLOVER_TYPES(Impl); + + template using iImplClover = iScalar, Ns>>; + typedef CompactWilsonCloverHelpers CompactHelpers; + + static void plusIdentity(const CloverField& in) { + int DimRep = Impl::Dimension; + + autoView(in_v, in, AcceleratorWrite); + + accelerator_for(ss, in.Grid()->oSites(), 1, { + for (int sa=0; saprec) { + NMAX++; + cond*=R/(double)(NMAX+1); + } + return NMAX; + } + + static int getNMAX(Lattice> &t, RealD R) {return getNMAX(1e-12,R);} + static int getNMAX(Lattice> &t, RealD R) {return getNMAX(1e-6,R);} + + static void MassTerm(CloverField& Clover, RealD diag_mass) { + // csw/(diag_mass) * clover + Clover = Clover * (1.0/diag_mass); + } + + static void Instantiate(CloverDiagonalField& Diagonal, CloverTriangleField& Triangle, + CloverDiagonalField& DiagonalInv, CloverTriangleField& TriangleInv, + RealD csw_t, RealD diag_mass) { + GridBase* grid = Diagonal.Grid(); + int NMAX = getNMAX(Diagonal, 3.*csw_t/diag_mass); + // To be optimized: too much memory traffic; implement exp in improved layout + // Felix + Fabian: replace code below with + // + // ConvertLayout Clover -> Diagonal, Triangle + // ModifyBoundaries + // EvaluateExp + + // code to be replaced + CloverField Clover(grid), ExpClover(grid); + + CompactHelpers::ConvertLayout(Diagonal, Triangle, Clover); + + ExpClover = Clover; + plusIdentity(ExpClover); // 1 + Clover + + RealD pref = 1.0; + for (int i=2; i<=NMAX; i++) { + Clover = Clover * Clover; + pref /= (RealD)(i); + ExpClover += pref * Clover; + } + + // Convert the data layout of the clover terms + CompactHelpers::ConvertLayout(ExpClover, Diagonal, Triangle); + Diagonal = Diagonal * diag_mass; + Triangle = Triangle * diag_mass; + CompactHelpers::ConvertLayout(adj(ExpClover), DiagonalInv, TriangleInv); + DiagonalInv = DiagonalInv*(1.0/diag_mass); + TriangleInv = TriangleInv*(1.0/diag_mass); + } + + + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) { + assert(0); + } + +}; + + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h index 3a166134..249b20bd 100644 --- a/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/CompactWilsonCloverFermion.h @@ -31,6 +31,7 @@ #include #include +#include NAMESPACE_BEGIN(Grid); @@ -85,7 +86,7 @@ NAMESPACE_BEGIN(Grid); // + (2 * 1 + 4 * 1/2) triangle parts = 4 triangle parts = 60 complex words per site // = 84 complex words per site -template +template class CompactWilsonCloverFermion : public WilsonFermion, public WilsonCloverHelpers, public CompactWilsonCloverHelpers { diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index 78cf7851..bfbc16b8 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -138,38 +138,62 @@ typedef WilsonTMFermion WilsonTMFermionF; typedef WilsonTMFermion WilsonTMFermionD; // Clover fermions -typedef WilsonCloverFermion WilsonCloverFermionR; -typedef WilsonCloverFermion WilsonCloverFermionF; -typedef WilsonCloverFermion WilsonCloverFermionD; +typedef CloverHelpers CloverR; +typedef CloverHelpers CloverF; +typedef CloverHelpers CloverD; -typedef WilsonCloverFermion WilsonCloverAdjFermionR; -typedef WilsonCloverFermion WilsonCloverAdjFermionF; -typedef WilsonCloverFermion WilsonCloverAdjFermionD; +typedef ExpCloverHelpers ExpCloverR; +typedef ExpCloverHelpers ExpCloverF; +typedef ExpCloverHelpers ExpCloverD; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionR; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionF; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionD; +typedef WilsonCloverFermion WilsonCloverFermionR; +typedef WilsonCloverFermion WilsonCloverFermionF; +typedef WilsonCloverFermion WilsonCloverFermionD; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionR; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionF; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionD; +typedef WilsonCloverFermion WilsonExpCloverFermionR; +typedef WilsonCloverFermion WilsonExpCloverFermionF; +typedef WilsonCloverFermion WilsonExpCloverFermionD; + +typedef WilsonCloverFermion WilsonCloverAdjFermionR; +typedef WilsonCloverFermion WilsonCloverAdjFermionF; +typedef WilsonCloverFermion WilsonCloverAdjFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionD; + +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionR; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionF; +typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionD; // Compact Clover fermions -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionD; +typedef CompactCloverHelpers CompactCloverR; +typedef CompactCloverHelpers CompactCloverF; +typedef CompactCloverHelpers CompactCloverD; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionD; +typedef CompactExpCloverHelpers CompactExpCloverR; +typedef CompactExpCloverHelpers CompactExpCloverF; +typedef CompactExpCloverHelpers CompactExpCloverD; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionD; +typedef CompactWilsonCloverFermion CompactWilsonCloverFermionR; +typedef CompactWilsonCloverFermion CompactWilsonCloverFermionF; +typedef CompactWilsonCloverFermion CompactWilsonCloverFermionD; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionD; +typedef CompactWilsonCloverFermion CompactWilsonExpCloverFermionR; +typedef CompactWilsonCloverFermion CompactWilsonExpCloverFermionF; +typedef CompactWilsonCloverFermion CompactWilsonExpCloverFermionD; + +typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionR; +typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionF; +typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionD; + +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionR; +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionF; +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionD; + +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionR; +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionF; +typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionD; // Domain Wall fermions typedef DomainWallFermion DomainWallFermionR; diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h index cd19fc10..562f08f4 100644 --- a/Grid/qcd/action/fermion/WilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h @@ -32,6 +32,7 @@ #include #include +#include NAMESPACE_BEGIN(Grid); @@ -51,7 +52,7 @@ NAMESPACE_BEGIN(Grid); // csw_r = csw_t to recover the isotropic version ////////////////////////////////////////////////////////////////// -template +template class WilsonCloverFermion : public WilsonFermion, public WilsonCloverHelpers { diff --git a/Grid/qcd/action/fermion/WilsonCloverHelpers.h b/Grid/qcd/action/fermion/WilsonCloverHelpers.h index 60f19317..c221d5f0 100644 --- a/Grid/qcd/action/fermion/WilsonCloverHelpers.h +++ b/Grid/qcd/action/fermion/WilsonCloverHelpers.h @@ -209,6 +209,8 @@ public: }; +//////////////////////////////////////////////////////// + template class CompactWilsonCloverHelpers { public: diff --git a/Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h index 3dfcb133..b42957ac 100644 --- a/Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/CompactWilsonCloverFermionImplementation.h @@ -32,17 +32,18 @@ #include #include + NAMESPACE_BEGIN(Grid); -template -CompactWilsonCloverFermion::CompactWilsonCloverFermion(GaugeField& _Umu, - GridCartesian& Fgrid, - GridRedBlackCartesian& Hgrid, - const RealD _mass, - const RealD _csw_r, - const RealD _csw_t, - const RealD _cF, - const WilsonAnisotropyCoefficients& clover_anisotropy, - const ImplParams& impl_p) +template +CompactWilsonCloverFermion::CompactWilsonCloverFermion(GaugeField& _Umu, + GridCartesian& Fgrid, + GridRedBlackCartesian& Hgrid, + const RealD _mass, + const RealD _csw_r, + const RealD _csw_t, + const RealD _cF, + const WilsonAnisotropyCoefficients& clover_anisotropy, + const ImplParams& impl_p) : WilsonBase(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy) , csw_r(_csw_r) , csw_t(_csw_t) @@ -68,40 +69,40 @@ CompactWilsonCloverFermion::CompactWilsonCloverFermion(GaugeField& _Umu, CompactHelpers::SetupMasks(this->BoundaryMask, this->BoundaryMaskEven, this->BoundaryMaskOdd); } -template -void CompactWilsonCloverFermion::Dhop(const FermionField& in, FermionField& out, int dag) { +template +void CompactWilsonCloverFermion::Dhop(const FermionField& in, FermionField& out, int dag) { WilsonBase::Dhop(in, out, dag); if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::DhopOE(const FermionField& in, FermionField& out, int dag) { +template +void CompactWilsonCloverFermion::DhopOE(const FermionField& in, FermionField& out, int dag) { WilsonBase::DhopOE(in, out, dag); if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::DhopEO(const FermionField& in, FermionField& out, int dag) { +template +void CompactWilsonCloverFermion::DhopEO(const FermionField& in, FermionField& out, int dag) { WilsonBase::DhopEO(in, out, dag); if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) { +template +void CompactWilsonCloverFermion::DhopDir(const FermionField& in, FermionField& out, int dir, int disp) { WilsonBase::DhopDir(in, out, dir, disp); if(this->open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::DhopDirAll(const FermionField& in, std::vector& out) { +template +void CompactWilsonCloverFermion::DhopDirAll(const FermionField& in, std::vector& out) { WilsonBase::DhopDirAll(in, out); if(this->open_boundaries) { for(auto& o : out) ApplyBoundaryMask(o); } } -template -void CompactWilsonCloverFermion::M(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::M(const FermionField& in, FermionField& out) { out.Checkerboard() = in.Checkerboard(); WilsonBase::Dhop(in, out, DaggerNo); // call base to save applying bc Mooee(in, Tmp); @@ -109,8 +110,8 @@ void CompactWilsonCloverFermion::M(const FermionField& in, FermionField& o if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::Mdag(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::Mdag(const FermionField& in, FermionField& out) { out.Checkerboard() = in.Checkerboard(); WilsonBase::Dhop(in, out, DaggerYes); // call base to save applying bc MooeeDag(in, Tmp); @@ -118,20 +119,20 @@ void CompactWilsonCloverFermion::Mdag(const FermionField& in, FermionField if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::Meooe(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::Meooe(const FermionField& in, FermionField& out) { WilsonBase::Meooe(in, out); if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::MeooeDag(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::MeooeDag(const FermionField& in, FermionField& out) { WilsonBase::MeooeDag(in, out); if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::Mooee(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::Mooee(const FermionField& in, FermionField& out) { if(in.Grid()->_isCheckerBoarded) { if(in.Checkerboard() == Odd) { MooeeInternal(in, out, DiagonalOdd, TriangleOdd); @@ -144,13 +145,13 @@ void CompactWilsonCloverFermion::Mooee(const FermionField& in, FermionFiel if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::MooeeDag(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::MooeeDag(const FermionField& in, FermionField& out) { Mooee(in, out); // blocks are hermitian } -template -void CompactWilsonCloverFermion::MooeeInv(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::MooeeInv(const FermionField& in, FermionField& out) { if(in.Grid()->_isCheckerBoarded) { if(in.Checkerboard() == Odd) { MooeeInternal(in, out, DiagonalInvOdd, TriangleInvOdd); @@ -163,23 +164,23 @@ void CompactWilsonCloverFermion::MooeeInv(const FermionField& in, FermionF if(open_boundaries) ApplyBoundaryMask(out); } -template -void CompactWilsonCloverFermion::MooeeInvDag(const FermionField& in, FermionField& out) { +template +void CompactWilsonCloverFermion::MooeeInvDag(const FermionField& in, FermionField& out) { MooeeInv(in, out); // blocks are hermitian } -template -void CompactWilsonCloverFermion::Mdir(const FermionField& in, FermionField& out, int dir, int disp) { +template +void CompactWilsonCloverFermion::Mdir(const FermionField& in, FermionField& out, int dir, int disp) { DhopDir(in, out, dir, disp); } -template -void CompactWilsonCloverFermion::MdirAll(const FermionField& in, std::vector& out) { +template +void CompactWilsonCloverFermion::MdirAll(const FermionField& in, std::vector& out) { DhopDirAll(in, out); } -template -void CompactWilsonCloverFermion::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) { +template +void CompactWilsonCloverFermion::MDeriv(GaugeField& force, const FermionField& X, const FermionField& Y, int dag) { assert(!open_boundaries); // TODO check for changes required for open bc // NOTE: code copied from original clover term @@ -251,7 +252,7 @@ void CompactWilsonCloverFermion::MDeriv(GaugeField& force, const FermionFi } PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok - force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu); // checked + force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked count++; } @@ -261,18 +262,18 @@ void CompactWilsonCloverFermion::MDeriv(GaugeField& force, const FermionFi force += clover_force; } -template -void CompactWilsonCloverFermion::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) { +template +void CompactWilsonCloverFermion::MooDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) { assert(0); } -template -void CompactWilsonCloverFermion::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) { +template +void CompactWilsonCloverFermion::MeeDeriv(GaugeField& mat, const FermionField& U, const FermionField& V, int dag) { assert(0); } -template -void CompactWilsonCloverFermion::MooeeInternal(const FermionField& in, +template +void CompactWilsonCloverFermion::MooeeInternal(const FermionField& in, FermionField& out, const CloverDiagonalField& diagonal, const CloverTriangleField& triangle) { @@ -285,8 +286,8 @@ void CompactWilsonCloverFermion::MooeeInternal(const FermionField& CompactHelpers::MooeeKernel(diagonal.oSites(), 1, in, out, diagonal, triangle); } -template -void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { +template +void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { // NOTE: parts copied from original implementation // Import gauge into base class @@ -318,8 +319,9 @@ void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { TmpOriginal += Helpers::fillCloverXT(Ex) * csw_t; TmpOriginal += Helpers::fillCloverYT(Ey) * csw_t; TmpOriginal += Helpers::fillCloverZT(Ez) * csw_t; - TmpOriginal += this->diag_mass; - + // Handle mass term based on clover policy + CloverHelpers::MassTerm(TmpOriginal, this->diag_mass); + // Convert the data layout of the clover term double t4 = usecond(); CompactHelpers::ConvertLayout(TmpOriginal, Diagonal, Triangle); @@ -328,9 +330,9 @@ void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { double t5 = usecond(); if(open_boundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, this->diag_mass); - // Invert the clover term in the improved layout + // Instantiate based on clover policy double t6 = usecond(); - CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv); + CloverHelpers::Instantiate(Diagonal, Triangle, DiagonalInv, TriangleInv, csw_t, this->diag_mass); // Fill the remaining clover fields double t7 = usecond(); @@ -353,7 +355,7 @@ void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { << ", fill clover = " << (t4 - t3) / 1e6 << ", convert = " << (t5 - t4) / 1e6 << ", boundaries = " << (t6 - t5) / 1e6 - << ", inversions = " << (t7 - t6) / 1e6 + << ", instantiate = " << (t7 - t6) / 1e6 << ", pick cbs = " << (t8 - t7) / 1e6 << ", total = " << (t8 - t0) / 1e6 << std::endl; diff --git a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h index e4d2e736..2ebf45f8 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h @@ -34,8 +34,8 @@ NAMESPACE_BEGIN(Grid); -template -WilsonCloverFermion::WilsonCloverFermion(GaugeField& _Umu, +template +WilsonCloverFermion::WilsonCloverFermion(GaugeField& _Umu, GridCartesian& Fgrid, GridRedBlackCartesian& Hgrid, const RealD _mass, @@ -74,8 +74,8 @@ WilsonCloverFermion::WilsonCloverFermion(GaugeField& } // *NOT* EO -template -void WilsonCloverFermion::M(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::M(const FermionField &in, FermionField &out) { FermionField temp(out.Grid()); @@ -89,8 +89,8 @@ void WilsonCloverFermion::M(const FermionField &in, FermionField &out) out += temp; } -template -void WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) { FermionField temp(out.Grid()); @@ -104,8 +104,8 @@ void WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) out += temp; } -template -void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) +template +void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) { double t0 = usecond(); WilsonFermion::ImportGauge(_Umu); @@ -131,47 +131,50 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) CloverTerm += Helpers::fillCloverXT(Ex) * csw_t; CloverTerm += Helpers::fillCloverYT(Ey) * csw_t; CloverTerm += Helpers::fillCloverZT(Ez) * csw_t; - CloverTerm += diag_mass; - + double t4 = usecond(); - int lvol = _Umu.Grid()->lSites(); - int DimRep = Impl::Dimension; + CloverHelpers::Instantiate(CloverTerm, CloverTermInv, csw_t, this->diag_mass); +// CloverTerm += diag_mass; +// +// double t4 = usecond(); +// int lvol = _Umu.Grid()->lSites(); +// int DimRep = Impl::Dimension; +// +// double t5 = usecond(); +// { +// autoView(CTv,CloverTerm,CpuRead); +// autoView(CTIv,CloverTermInv,CpuWrite); +// thread_for(site, lvol, { +// Coordinate lcoor; +// grid->LocalIndexToLocalCoor(site, lcoor); +// Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); +// Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); +// typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero(); +// peekLocalSite(Qx, CTv, lcoor); +// //if (csw!=0){ +// for (int j = 0; j < Ns; j++) +// for (int k = 0; k < Ns; k++) +// for (int a = 0; a < DimRep; a++) +// for (int b = 0; b < DimRep; b++){ +// auto zz = Qx()(j, k)(a, b); +// EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex(zz); +// } +// // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; +// +// EigenInvCloverOp = EigenCloverOp.inverse(); +// //std::cout << EigenInvCloverOp << std::endl; +// for (int j = 0; j < Ns; j++) +// for (int k = 0; k < Ns; k++) +// for (int a = 0; a < DimRep; a++) +// for (int b = 0; b < DimRep; b++) +// Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); +// // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; +// // } +// pokeLocalSite(Qxinv, CTIv, lcoor); +// }); +// } double t5 = usecond(); - { - autoView(CTv,CloverTerm,CpuRead); - autoView(CTIv,CloverTermInv,CpuWrite); - thread_for(site, lvol, { - Coordinate lcoor; - grid->LocalIndexToLocalCoor(site, lcoor); - Eigen::MatrixXcd EigenCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); - Eigen::MatrixXcd EigenInvCloverOp = Eigen::MatrixXcd::Zero(Ns * DimRep, Ns * DimRep); - typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero(); - peekLocalSite(Qx, CTv, lcoor); - //if (csw!=0){ - for (int j = 0; j < Ns; j++) - for (int k = 0; k < Ns; k++) - for (int a = 0; a < DimRep; a++) - for (int b = 0; b < DimRep; b++){ - auto zz = Qx()(j, k)(a, b); - EigenCloverOp(a + j * DimRep, b + k * DimRep) = std::complex(zz); - } - // if (site==0) std::cout << "site =" << site << "\n" << EigenCloverOp << std::endl; - - EigenInvCloverOp = EigenCloverOp.inverse(); - //std::cout << EigenInvCloverOp << std::endl; - for (int j = 0; j < Ns; j++) - for (int k = 0; k < Ns; k++) - for (int a = 0; a < DimRep; a++) - for (int b = 0; b < DimRep; b++) - Qxinv()(j, k)(a, b) = EigenInvCloverOp(a + j * DimRep, b + k * DimRep); - // if (site==0) std::cout << "site =" << site << "\n" << EigenInvCloverOp << std::endl; - // } - pokeLocalSite(Qxinv, CTIv, lcoor); - }); - } - - double t6 = usecond(); // Separate the even and odd parts pickCheckerboard(Even, CloverTermEven, CloverTerm); pickCheckerboard(Odd, CloverTermOdd, CloverTerm); @@ -184,7 +187,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); - double t7 = usecond(); + double t6 = usecond(); #if 0 std::cout << GridLogMessage << "WilsonCloverFermion::ImportGauge timings:" @@ -192,40 +195,39 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) << ", allocations = " << (t2 - t1) / 1e6 << ", field strength = " << (t3 - t2) / 1e6 << ", fill clover = " << (t4 - t3) / 1e6 - << ", misc = " << (t5 - t4) / 1e6 - << ", inversions = " << (t6 - t5) / 1e6 - << ", pick cbs = " << (t7 - t6) / 1e6 - << ", total = " << (t7 - t0) / 1e6 + << ", instantiation = " << (t5 - t4) / 1e6 + << ", pick cbs = " << (t6 - t5) / 1e6 + << ", total = " << (t6 - t0) / 1e6 << std::endl; #endif } -template -void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::Mooee(const FermionField &in, FermionField &out) { this->MooeeInternal(in, out, DaggerNo, InverseNo); } -template -void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::MooeeDag(const FermionField &in, FermionField &out) { this->MooeeInternal(in, out, DaggerYes, InverseNo); } -template -void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::MooeeInv(const FermionField &in, FermionField &out) { this->MooeeInternal(in, out, DaggerNo, InverseYes); } -template -void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) +template +void WilsonCloverFermion::MooeeInvDag(const FermionField &in, FermionField &out) { this->MooeeInternal(in, out, DaggerYes, InverseYes); } -template -void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv) +template +void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv) { out.Checkerboard() = in.Checkerboard(); CloverField *Clover; @@ -278,8 +280,8 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie } // MooeeInternal // Derivative parts unpreconditioned pseudofermions -template -void WilsonCloverFermion::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) +template +void WilsonCloverFermion::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) { conformable(X.Grid(), Y.Grid()); conformable(X.Grid(), force.Grid()); @@ -349,7 +351,7 @@ void WilsonCloverFermion::MDeriv(GaugeField &force, const FermionField &X, } PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok - force_mu -= factor*Helpers::Cmunu(U, lambda, mu, nu); // checked + force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked count++; } @@ -360,15 +362,15 @@ void WilsonCloverFermion::MDeriv(GaugeField &force, const FermionField &X, } // Derivative parts -template -void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) +template +void WilsonCloverFermion::MooDeriv(GaugeField &mat, const FermionField &X, const FermionField &Y, int dag) { assert(0); } // Derivative parts -template -void WilsonCloverFermion::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) +template +void WilsonCloverFermion::MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) { assert(0); // not implemented yet } diff --git a/Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master b/Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master index 7c91c252..1ea1549c 100644 --- a/Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master +++ b/Grid/qcd/action/fermion/instantiation/CompactWilsonCloverFermionInstantiation.cc.master @@ -9,6 +9,7 @@ Author: paboyle Author: Guido Cossu Author: Daniel Richtmann + Author: Mattia Bruno 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 @@ -32,10 +33,12 @@ #include #include #include +#include NAMESPACE_BEGIN(Grid); #include "impl.h" -template class CompactWilsonCloverFermion; +template class CompactWilsonCloverFermion>; +template class CompactWilsonCloverFermion>; NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master b/Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master index af99dfb6..2cea3656 100644 --- a/Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master +++ b/Grid/qcd/action/fermion/instantiation/WilsonCloverFermionInstantiation.cc.master @@ -8,7 +8,8 @@ Author: paboyle Author: Guido Cossu - + Author: Mattia Bruno + 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 @@ -31,10 +32,12 @@ #include #include #include +#include NAMESPACE_BEGIN(Grid); #include "impl.h" -template class WilsonCloverFermion; +template class WilsonCloverFermion>; +template class WilsonCloverFermion>; NAMESPACE_END(Grid);