diff --git a/Grid/qcd/action/fermion/WilsonCloverFermion.h b/Grid/qcd/action/fermion/WilsonCloverFermion.h index 92af7111..cd19fc10 100644 --- a/Grid/qcd/action/fermion/WilsonCloverFermion.h +++ b/Grid/qcd/action/fermion/WilsonCloverFermion.h @@ -4,10 +4,11 @@ Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.h - Copyright (C) 2017 + Copyright (C) 2017 - 2022 Author: Guido Cossu Author: David Preti <> + Author: Daniel Richtmann 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 @@ -29,7 +30,8 @@ #pragma once -#include +#include +#include NAMESPACE_BEGIN(Grid); @@ -50,18 +52,15 @@ NAMESPACE_BEGIN(Grid); ////////////////////////////////////////////////////////////////// template -class WilsonCloverFermion : public WilsonFermion +class WilsonCloverFermion : public WilsonFermion, + public WilsonCloverHelpers { public: - // Types definitions INHERIT_IMPL_TYPES(Impl); - template - using iImplClover = iScalar, Ns>>; - typedef iImplClover SiteCloverType; - typedef Lattice CloverFieldType; + INHERIT_CLOVER_TYPES(Impl); -public: - typedef WilsonFermion WilsonBase; + typedef WilsonFermion WilsonBase; + typedef WilsonCloverHelpers Helpers; virtual int ConstEE(void) { return 0; }; virtual void Instantiatable(void){}; @@ -72,42 +71,7 @@ public: const RealD _csw_r = 0.0, const RealD _csw_t = 0.0, const WilsonAnisotropyCoefficients &clover_anisotropy = WilsonAnisotropyCoefficients(), - const ImplParams &impl_p = ImplParams()) : WilsonFermion(_Umu, - Fgrid, - Hgrid, - _mass, impl_p, clover_anisotropy), - CloverTerm(&Fgrid), - CloverTermInv(&Fgrid), - CloverTermEven(&Hgrid), - CloverTermOdd(&Hgrid), - CloverTermInvEven(&Hgrid), - CloverTermInvOdd(&Hgrid), - CloverTermDagEven(&Hgrid), - CloverTermDagOdd(&Hgrid), - CloverTermInvDagEven(&Hgrid), - CloverTermInvDagOdd(&Hgrid) - { - assert(Nd == 4); // require 4 dimensions - - if (clover_anisotropy.isAnisotropic) - { - csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0; - diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); - } - else - { - csw_r = _csw_r * 0.5; - diag_mass = 4.0 + _mass; - } - csw_t = _csw_t * 0.5; - - if (csw_r == 0) - std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl; - if (csw_t == 0) - std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl; - - ImportGauge(_Umu); - } + const ImplParams &impl_p = ImplParams()); virtual void M(const FermionField &in, FermionField &out); virtual void Mdag(const FermionField &in, FermionField &out); @@ -124,250 +88,21 @@ public: void ImportGauge(const GaugeField &_Umu); // Derivative parts unpreconditioned pseudofermions - void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) - { - conformable(X.Grid(), Y.Grid()); - conformable(X.Grid(), force.Grid()); - GaugeLinkField force_mu(force.Grid()), lambda(force.Grid()); - GaugeField clover_force(force.Grid()); - PropagatorField Lambda(force.Grid()); + void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag); - // Guido: Here we are hitting some performance issues: - // need to extract the components of the DoubledGaugeField - // for each call - // Possible solution - // Create a vector object to store them? (cons: wasting space) - std::vector U(Nd, this->Umu.Grid()); - - Impl::extractLinkField(U, this->Umu); - - force = Zero(); - // Derivative of the Wilson hopping term - this->DhopDeriv(force, X, Y, dag); - - /////////////////////////////////////////////////////////// - // 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*Cmunu(U, lambda, mu, nu); // checked - count++; - } - - pokeLorentz(clover_force, U[mu] * force_mu, mu); - } - //clover_force *= csw; - force += clover_force; - } - - // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis - GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) - { - conformable(lambda.Grid(), U[0].Grid()); - GaugeLinkField out(lambda.Grid()), tmp(lambda.Grid()); - // insertion in upper staple - // please check redundancy of shift operations - - // C1+ - tmp = lambda * U[nu]; - out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - - // C2+ - tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu); - out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - - // C3+ - tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu); - out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); - - // C4+ - out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda; - - // insertion in lower staple - // C1- - out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); - - // C2- - tmp = adj(lambda) * U[nu]; - out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); - - // C3- - tmp = lambda * U[nu]; - out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); - - // C4- - out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda; - - return out; - } - -protected: +public: // here fixing the 4 dimensions, make it more general? RealD csw_r; // Clover coefficient - spatial RealD csw_t; // Clover coefficient - temporal RealD diag_mass; // Mass term - CloverFieldType CloverTerm, CloverTermInv; // Clover term - CloverFieldType CloverTermEven, CloverTermOdd; // Clover term EO - CloverFieldType CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO - CloverFieldType CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO - CloverFieldType CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO - - public: - // eventually these can be compressed into 6x6 blocks instead of the 12x12 - // using the DeGrand-Rossi basis for the gamma matrices - CloverFieldType fillCloverYZ(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - T = Zero(); - autoView(T_v,T,AcceleratorWrite); - autoView(F_v,F,AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 1) = timesMinusI(F_v[i]()()); - T_v[i]()(1, 0) = timesMinusI(F_v[i]()()); - T_v[i]()(2, 3) = timesMinusI(F_v[i]()()); - T_v[i]()(3, 2) = timesMinusI(F_v[i]()()); - }); - - return T; - } - - CloverFieldType fillCloverXZ(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - T = Zero(); - - autoView(T_v, T,AcceleratorWrite); - autoView(F_v, F,AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 1) = -F_v[i]()(); - T_v[i]()(1, 0) = F_v[i]()(); - T_v[i]()(2, 3) = -F_v[i]()(); - T_v[i]()(3, 2) = F_v[i]()(); - }); - - return T; - } - - CloverFieldType fillCloverXY(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - T = Zero(); - - autoView(T_v,T,AcceleratorWrite); - autoView(F_v,F,AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 0) = timesMinusI(F_v[i]()()); - T_v[i]()(1, 1) = timesI(F_v[i]()()); - T_v[i]()(2, 2) = timesMinusI(F_v[i]()()); - T_v[i]()(3, 3) = timesI(F_v[i]()()); - }); - - return T; - } - - CloverFieldType fillCloverXT(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - T = Zero(); - - autoView( T_v , T, AcceleratorWrite); - autoView( F_v , F, AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 1) = timesI(F_v[i]()()); - T_v[i]()(1, 0) = timesI(F_v[i]()()); - T_v[i]()(2, 3) = timesMinusI(F_v[i]()()); - T_v[i]()(3, 2) = timesMinusI(F_v[i]()()); - }); - - return T; - } - - CloverFieldType fillCloverYT(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - T = Zero(); - - autoView( T_v ,T,AcceleratorWrite); - autoView( F_v ,F,AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 1) = -(F_v[i]()()); - T_v[i]()(1, 0) = (F_v[i]()()); - T_v[i]()(2, 3) = (F_v[i]()()); - T_v[i]()(3, 2) = -(F_v[i]()()); - }); - - return T; - } - - CloverFieldType fillCloverZT(const GaugeLinkField &F) - { - CloverFieldType T(F.Grid()); - - T = Zero(); - - autoView( T_v , T,AcceleratorWrite); - autoView( F_v , F,AcceleratorRead); - accelerator_for(i, CloverTerm.Grid()->oSites(),1, - { - T_v[i]()(0, 0) = timesI(F_v[i]()()); - T_v[i]()(1, 1) = timesMinusI(F_v[i]()()); - T_v[i]()(2, 2) = timesMinusI(F_v[i]()()); - T_v[i]()(3, 3) = timesI(F_v[i]()()); - }); - - return T; - } + CloverField CloverTerm, CloverTermInv; // Clover term + CloverField CloverTermEven, CloverTermOdd; // Clover term EO + CloverField CloverTermInvEven, CloverTermInvOdd; // Clover term Inv EO + CloverField CloverTermDagEven, CloverTermDagOdd; // Clover term Dag EO + CloverField CloverTermInvDagEven, CloverTermInvDagOdd; // Clover term Inv Dag EO }; + NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/WilsonCloverHelpers.h b/Grid/qcd/action/fermion/WilsonCloverHelpers.h new file mode 100644 index 00000000..80f9032c --- /dev/null +++ b/Grid/qcd/action/fermion/WilsonCloverHelpers.h @@ -0,0 +1,191 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverHelpers.h + + Copyright (C) 2021 - 2022 + + Author: Daniel Richtmann + + 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 + +// Helper routines that implement common clover functionality + +NAMESPACE_BEGIN(Grid); + +template class WilsonCloverHelpers { +public: + + INHERIT_IMPL_TYPES(Impl); + INHERIT_CLOVER_TYPES(Impl); + + // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) + { + conformable(lambda.Grid(), U[0].Grid()); + GaugeLinkField out(lambda.Grid()), tmp(lambda.Grid()); + // insertion in upper staple + // please check redundancy of shift operations + + // C1+ + tmp = lambda * U[nu]; + out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C2+ + tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); + + // C3+ + tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu); + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); + + // C4+ + out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda; + + // insertion in lower staple + // C1- + out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C2- + tmp = adj(lambda) * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu); + + // C3- + tmp = lambda * U[nu]; + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu); + + // C4- + out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda; + + return out; + } + + static CloverField fillCloverYZ(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + T = Zero(); + autoView(T_v,T,AcceleratorWrite); + autoView(F_v,F,AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 1) = timesMinusI(F_v[i]()()); + T_v[i]()(1, 0) = timesMinusI(F_v[i]()()); + T_v[i]()(2, 3) = timesMinusI(F_v[i]()()); + T_v[i]()(3, 2) = timesMinusI(F_v[i]()()); + }); + + return T; + } + + static CloverField fillCloverXZ(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + T = Zero(); + + autoView(T_v, T,AcceleratorWrite); + autoView(F_v, F,AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 1) = -F_v[i]()(); + T_v[i]()(1, 0) = F_v[i]()(); + T_v[i]()(2, 3) = -F_v[i]()(); + T_v[i]()(3, 2) = F_v[i]()(); + }); + + return T; + } + + static CloverField fillCloverXY(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + T = Zero(); + + autoView(T_v,T,AcceleratorWrite); + autoView(F_v,F,AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 0) = timesMinusI(F_v[i]()()); + T_v[i]()(1, 1) = timesI(F_v[i]()()); + T_v[i]()(2, 2) = timesMinusI(F_v[i]()()); + T_v[i]()(3, 3) = timesI(F_v[i]()()); + }); + + return T; + } + + static CloverField fillCloverXT(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + T = Zero(); + + autoView( T_v , T, AcceleratorWrite); + autoView( F_v , F, AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 1) = timesI(F_v[i]()()); + T_v[i]()(1, 0) = timesI(F_v[i]()()); + T_v[i]()(2, 3) = timesMinusI(F_v[i]()()); + T_v[i]()(3, 2) = timesMinusI(F_v[i]()()); + }); + + return T; + } + + static CloverField fillCloverYT(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + T = Zero(); + + autoView( T_v ,T,AcceleratorWrite); + autoView( F_v ,F,AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 1) = -(F_v[i]()()); + T_v[i]()(1, 0) = (F_v[i]()()); + T_v[i]()(2, 3) = (F_v[i]()()); + T_v[i]()(3, 2) = -(F_v[i]()()); + }); + + return T; + } + + static CloverField fillCloverZT(const GaugeLinkField &F) + { + CloverField T(F.Grid()); + + T = Zero(); + + autoView( T_v , T,AcceleratorWrite); + autoView( F_v , F,AcceleratorRead); + accelerator_for(i, T.Grid()->oSites(),1, + { + T_v[i]()(0, 0) = timesI(F_v[i]()()); + T_v[i]()(1, 1) = timesMinusI(F_v[i]()()); + T_v[i]()(2, 2) = timesMinusI(F_v[i]()()); + T_v[i]()(3, 3) = timesI(F_v[i]()()); + }); + + return T; + } +}; + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/WilsonCloverTypes.h b/Grid/qcd/action/fermion/WilsonCloverTypes.h new file mode 100644 index 00000000..4428259f --- /dev/null +++ b/Grid/qcd/action/fermion/WilsonCloverTypes.h @@ -0,0 +1,49 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/action/fermion/WilsonCloverTypes.h + + Copyright (C) 2021 - 2022 + + Author: Daniel Richtmann + + 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); + +template +class WilsonCloverTypes { +public: + INHERIT_IMPL_TYPES(Impl); + + template using iImplClover = iScalar, Ns>>; + + typedef iImplClover SiteClover; + + typedef Lattice CloverField; +}; + +#define INHERIT_CLOVER_TYPES(Impl) \ + typedef typename WilsonCloverTypes::SiteClover SiteClover; \ + typedef typename WilsonCloverTypes::CloverField CloverField; + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h index 3032a80c..fd81af11 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h @@ -2,12 +2,13 @@ Grid physics library, www.github.com/paboyle/Grid - Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc + Source file: ./lib/qcd/action/fermion/WilsonCloverFermionImplementation.h - Copyright (C) 2017 + Copyright (C) 2017 - 2022 Author: paboyle Author: Guido Cossu + Author: Daniel Richtmann 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 @@ -33,6 +34,45 @@ NAMESPACE_BEGIN(Grid); +template +WilsonCloverFermion::WilsonCloverFermion(GaugeField& _Umu, + GridCartesian& Fgrid, + GridRedBlackCartesian& Hgrid, + const RealD _mass, + const RealD _csw_r, + const RealD _csw_t, + const WilsonAnisotropyCoefficients& clover_anisotropy, + const ImplParams& impl_p) + : WilsonFermion(_Umu, Fgrid, Hgrid, _mass, impl_p, clover_anisotropy) + , CloverTerm(&Fgrid) + , CloverTermInv(&Fgrid) + , CloverTermEven(&Hgrid) + , CloverTermOdd(&Hgrid) + , CloverTermInvEven(&Hgrid) + , CloverTermInvOdd(&Hgrid) + , CloverTermDagEven(&Hgrid) + , CloverTermDagOdd(&Hgrid) + , CloverTermInvDagEven(&Hgrid) + , CloverTermInvDagOdd(&Hgrid) { + assert(Nd == 4); // require 4 dimensions + + if(clover_anisotropy.isAnisotropic) { + csw_r = _csw_r * 0.5 / clover_anisotropy.xi_0; + diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); + } else { + csw_r = _csw_r * 0.5; + diag_mass = 4.0 + _mass; + } + csw_t = _csw_t * 0.5; + + if(csw_r == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_r = 0" << std::endl; + if(csw_t == 0) + std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw_t = 0" << std::endl; + + ImportGauge(_Umu); +} + // *NOT* EO template void WilsonCloverFermion::M(const FermionField &in, FermionField &out) @@ -67,10 +107,13 @@ void WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) template void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) { + double t0 = usecond(); WilsonFermion::ImportGauge(_Umu); + double t1 = usecond(); GridBase *grid = _Umu.Grid(); typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); + double t2 = usecond(); // Compute the field strength terms mu>nu WilsonLoops::FieldStrength(Bx, _Umu, Zdir, Ydir); WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); @@ -79,19 +122,22 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) WilsonLoops::FieldStrength(Ey, _Umu, Tdir, Ydir); WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); + double t3 = usecond(); // Compute the Clover Operator acting on Colour and Spin // multiply here by the clover coefficients for the anisotropy - CloverTerm = fillCloverYZ(Bx) * csw_r; - CloverTerm += fillCloverXZ(By) * csw_r; - CloverTerm += fillCloverXY(Bz) * csw_r; - CloverTerm += fillCloverXT(Ex) * csw_t; - CloverTerm += fillCloverYT(Ey) * csw_t; - CloverTerm += fillCloverZT(Ez) * csw_t; + CloverTerm = Helpers::fillCloverYZ(Bx) * csw_r; + CloverTerm += Helpers::fillCloverXZ(By) * csw_r; + CloverTerm += Helpers::fillCloverXY(Bz) * csw_r; + 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; + double t5 = usecond(); { autoView(CTv,CloverTerm,CpuRead); autoView(CTIv,CloverTermInv,CpuWrite); @@ -100,7 +146,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) 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 SiteCloverType::scalar_object Qx = Zero(), Qxinv = Zero(); + typename SiteClover::scalar_object Qx = Zero(), Qxinv = Zero(); peekLocalSite(Qx, CTv, lcoor); //if (csw!=0){ for (int j = 0; j < Ns; j++) @@ -125,6 +171,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) }); } + double t6 = usecond(); // Separate the even and odd parts pickCheckerboard(Even, CloverTermEven, CloverTerm); pickCheckerboard(Odd, CloverTermOdd, CloverTerm); @@ -137,6 +184,20 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) pickCheckerboard(Even, CloverTermInvDagEven, adj(CloverTermInv)); pickCheckerboard(Odd, CloverTermInvDagOdd, adj(CloverTermInv)); + double t7 = usecond(); + +#if 0 + std::cout << GridLogMessage << "WilsonCloverFermion::ImportGauge timings:" + << " WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 + << ", 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 + << std::endl; +#endif } template @@ -167,7 +228,7 @@ template void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv) { out.Checkerboard() = in.Checkerboard(); - CloverFieldType *Clover; + CloverField *Clover; assert(in.Checkerboard() == Odd || in.Checkerboard() == Even); if (dag) @@ -214,9 +275,89 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie out = *Clover * in; } } - } // MooeeInternal +// Derivative parts unpreconditioned pseudofermions +template +void WilsonCloverFermion::MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag) +{ + conformable(X.Grid(), Y.Grid()); + conformable(X.Grid(), force.Grid()); + GaugeLinkField force_mu(force.Grid()), lambda(force.Grid()); + GaugeField clover_force(force.Grid()); + PropagatorField Lambda(force.Grid()); + + // Guido: Here we are hitting some performance issues: + // need to extract the components of the DoubledGaugeField + // for each call + // Possible solution + // Create a vector object to store them? (cons: wasting space) + std::vector U(Nd, this->Umu.Grid()); + + Impl::extractLinkField(U, this->Umu); + + force = Zero(); + // Derivative of the Wilson hopping term + this->DhopDeriv(force, X, Y, dag); + + /////////////////////////////////////////////////////////// + // 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); + } + //clover_force *= csw; + force += clover_force; +} // Derivative parts template