diff --git a/Grid/DisableWarnings.h b/Grid/DisableWarnings.h index 4bd1edd0..64e4faf4 100644 --- a/Grid/DisableWarnings.h +++ b/Grid/DisableWarnings.h @@ -44,14 +44,22 @@ directory #ifdef __NVCC__ //disables nvcc specific warning in json.hpp #pragma clang diagnostic ignored "-Wdeprecated-register" + +#if (__CUDACC_VER_MAJOR__ >= 11) && (__CUDACC_VER_MINOR__ >= 5) + //disables nvcc specific warning in json.hpp +#pragma nv_diag_suppress unsigned_compare_with_zero +#pragma nv_diag_suppress cast_to_qualified_type + //disables nvcc specific warning in many files +#pragma nv_diag_suppress esa_on_defaulted_function_ignored +#pragma nv_diag_suppress extra_semicolon +#else + //disables nvcc specific warning in json.hpp #pragma diag_suppress unsigned_compare_with_zero #pragma diag_suppress cast_to_qualified_type - //disables nvcc specific warning in many files #pragma diag_suppress esa_on_defaulted_function_ignored #pragma diag_suppress extra_semicolon - -//Eigen only +#endif #endif // Disable vectorisation in Eigen on the Power8/9 and PowerPC diff --git a/Grid/Grid_Eigen_Dense.h b/Grid/Grid_Eigen_Dense.h index c62d9cdb..5aee81de 100644 --- a/Grid/Grid_Eigen_Dense.h +++ b/Grid/Grid_Eigen_Dense.h @@ -14,7 +14,11 @@ /* NVCC save and restore compile environment*/ #ifdef __NVCC__ #pragma push +#if (__CUDACC_VER_MAJOR__ >= 11) && (__CUDACC_VER_MINOR__ >= 5) +#pragma nv_diag_suppress code_is_unreachable +#else #pragma diag_suppress code_is_unreachable +#endif #pragma push_macro("__CUDA_ARCH__") #pragma push_macro("__NVCC__") #pragma push_macro("__CUDACC__") diff --git a/Grid/algorithms/iterative/Deflation.h b/Grid/algorithms/iterative/Deflation.h index 2eb28bf9..1a8f97c9 100644 --- a/Grid/algorithms/iterative/Deflation.h +++ b/Grid/algorithms/iterative/Deflation.h @@ -113,7 +113,43 @@ public: blockPromote(guess_coarse,guess,subspace); guess.Checkerboard() = src.Checkerboard(); }; -}; + + void operator()(const std::vector &src,std::vector &guess) { + int Nevec = (int)evec_coarse.size(); + int Nsrc = (int)src.size(); + // make temp variables + std::vector src_coarse(Nsrc,evec_coarse[0].Grid()); + std::vector guess_coarse(Nsrc,evec_coarse[0].Grid()); + //Preporcessing + std::cout << GridLogMessage << "Start BlockProject for loop" << std::endl; + for (int j=0;j #include #include +#include #include #include @@ -654,7 +655,8 @@ class IldgWriter : public ScidacWriter { // Fill ILDG header data struct ////////////////////////////////////////////////////// ildgFormat ildgfmt ; - ildgfmt.field = std::string("su3gauge"); + const std::string stNC = std::to_string( Nc ) ; + ildgfmt.field = std::string("su"+stNC+"gauge"); if ( format == std::string("IEEE32BIG") ) { ildgfmt.precision = 32; @@ -871,7 +873,8 @@ class IldgReader : public GridLimeReader { } else { assert(found_ildgFormat); - assert ( ildgFormat_.field == std::string("su3gauge") ); + const std::string stNC = std::to_string( Nc ) ; + assert ( ildgFormat_.field == std::string("su"+stNC+"gauge") ); /////////////////////////////////////////////////////////////////////////////////////// // Populate our Grid metadata as best we can @@ -879,7 +882,7 @@ class IldgReader : public GridLimeReader { std::ostringstream vers; vers << ildgFormat_.version; FieldMetaData_.hdr_version = vers.str(); - FieldMetaData_.data_type = std::string("4D_SU3_GAUGE_3X3"); + FieldMetaData_.data_type = std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC); FieldMetaData_.nd=4; FieldMetaData_.dimension.resize(4); diff --git a/Grid/parallelIO/MetaData.h b/Grid/parallelIO/MetaData.h index af8b3f76..6b9d8708 100644 --- a/Grid/parallelIO/MetaData.h +++ b/Grid/parallelIO/MetaData.h @@ -6,8 +6,8 @@ Copyright (C) 2015 - Author: Peter Boyle + Author: Jamie Hudspith 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 @@ -182,8 +182,8 @@ class GaugeStatistics public: void operator()(Lattice & data,FieldMetaData &header) { - header.link_trace=WilsonLoops::linkTrace(data); - header.plaquette =WilsonLoops::avgPlaquette(data); + header.link_trace = WilsonLoops::linkTrace(data); + header.plaquette = WilsonLoops::avgPlaquette(data); } }; typedef GaugeStatistics PeriodicGaugeStatistics; @@ -203,20 +203,24 @@ template<> inline void PrepareMetaData(Lattice 1 ) ; for(int mu=0;mu using iLorentzColour2x3 = iVector, 2>, Nd >; +template using iLorentzColour2x3 = iVector, Nc-1>, Nd >; typedef iLorentzColour2x3 LorentzColour2x3; typedef iLorentzColour2x3 LorentzColour2x3F; @@ -278,7 +282,6 @@ struct GaugeSimpleMunger{ template struct GaugeSimpleUnmunger { - void operator()(sobj &in, fobj &out) { for (int mu = 0; mu < Nd; mu++) { for (int i = 0; i < Nc; i++) { @@ -317,8 +320,8 @@ template struct Gauge3x2munger{ void operator() (fobj &in,sobj &out){ for(int mu=0;mu struct Gauge3x2unmunger{ void operator() (sobj &in,fobj &out){ for(int mu=0;mu Author: Peter Boyle Author: paboyle + Author: Jamie Hudspith 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 @@ -30,6 +31,8 @@ #ifndef GRID_NERSC_IO_H #define GRID_NERSC_IO_H +#include + NAMESPACE_BEGIN(Grid); using namespace Grid; @@ -145,15 +148,17 @@ public: std::string format(header.floating_point); - int ieee32big = (format == std::string("IEEE32BIG")); - int ieee32 = (format == std::string("IEEE32")); - int ieee64big = (format == std::string("IEEE64BIG")); - int ieee64 = (format == std::string("IEEE64") || format == std::string("IEEE64LITTLE")); + const int ieee32big = (format == std::string("IEEE32BIG")); + const int ieee32 = (format == std::string("IEEE32")); + const int ieee64big = (format == std::string("IEEE64BIG")); + const int ieee64 = (format == std::string("IEEE64") || \ + format == std::string("IEEE64LITTLE")); uint32_t nersc_csum,scidac_csuma,scidac_csumb; // depending on datatype, set up munger; // munger is a function of - if ( header.data_type == std::string("4D_SU3_GAUGE") ) { + const std::string stNC = std::to_string( Nc ) ; + if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE") ) { if ( ieee32 || ieee32big ) { BinaryIO::readLatticeObject (Umu,file,Gauge3x2munger(), offset,format, @@ -164,7 +169,7 @@ public: (Umu,file,Gauge3x2munger(),offset,format, nersc_csum,scidac_csuma,scidac_csumb); } - } else if ( header.data_type == std::string("4D_SU3_GAUGE_3x3") ) { + } else if ( header.data_type == std::string("4D_SU"+stNC+"_GAUGE_"+stNC+"x"+stNC) ) { if ( ieee32 || ieee32big ) { BinaryIO::readLatticeObject (Umu,file,GaugeSimpleMunger(),offset,format, @@ -209,27 +214,29 @@ public: template static inline void writeConfiguration(Lattice &Umu, std::string file, - std::string ens_label = std::string("DWF")) + std::string ens_label = std::string("DWF"), + std::string ens_id = std::string("UKQCD"), + unsigned int sequence_number = 1) { - writeConfiguration(Umu,file,0,1,ens_label); + writeConfiguration(Umu,file,0,1,ens_label,ens_id,sequence_number); } template static inline void writeConfiguration(Lattice &Umu, std::string file, int two_row, int bits32, - std::string ens_label = std::string("DWF")) + std::string ens_label = std::string("DWF"), + std::string ens_id = std::string("UKQCD"), + unsigned int sequence_number = 1) { typedef vLorentzColourMatrixD vobj; typedef typename vobj::scalar_object sobj; FieldMetaData header; - /////////////////////////////////////////// - // Following should become arguments - /////////////////////////////////////////// - header.sequence_number = 1; - header.ensemble_id = std::string("UKQCD"); + header.sequence_number = sequence_number; + header.ensemble_id = ens_id; header.ensemble_label = ens_label; + header.hdr_version = "1.0" ; typedef LorentzColourMatrixD fobj3D; typedef LorentzColour2x3D fobj2D; @@ -243,10 +250,14 @@ public: uint64_t offset; - // Sod it -- always write 3x3 double - header.floating_point = std::string("IEEE64BIG"); - header.data_type = std::string("4D_SU3_GAUGE_3x3"); - GaugeSimpleUnmunger munge; + // Sod it -- always write NcxNc double + header.floating_point = std::string("IEEE64BIG"); + const std::string stNC = std::to_string( Nc ) ; + if( two_row ) { + header.data_type = std::string("4D_SU" + stNC + "_GAUGE" ); + } else { + header.data_type = std::string("4D_SU" + stNC + "_GAUGE_" + stNC + "x" + stNC ); + } if ( grid->IsBoss() ) { truncate(file); offset = writeHeader(header,file); @@ -254,8 +265,15 @@ public: grid->Broadcast(0,(void *)&offset,sizeof(offset)); uint32_t nersc_csum,scidac_csuma,scidac_csumb; - BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, - nersc_csum,scidac_csuma,scidac_csumb); + if( two_row ) { + Gauge3x2unmunger munge; + BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, + nersc_csum,scidac_csuma,scidac_csumb); + } else { + GaugeSimpleUnmunger munge; + BinaryIO::writeLatticeObject(Umu,file,munge,offset,header.floating_point, + nersc_csum,scidac_csuma,scidac_csumb); + } header.checksum = nersc_csum; if ( grid->IsBoss() ) { writeHeader(header,file); @@ -287,8 +305,7 @@ public: header.plaquette=0.0; MachineCharacteristics(header); - uint64_t offset; - + uint64_t offset; #ifdef RNG_RANLUX header.floating_point = std::string("UINT64"); header.data_type = std::string("RANLUX48"); @@ -328,7 +345,7 @@ public: GridBase *grid = parallel.Grid(); - uint64_t offset = readHeader(file,grid,header); + uint64_t offset = readHeader(file,grid,header); FieldMetaData clone(header); diff --git a/Grid/pugixml/pugixml.cc b/Grid/pugixml/pugixml.cc index 45e6496a..81973964 100644 --- a/Grid/pugixml/pugixml.cc +++ b/Grid/pugixml/pugixml.cc @@ -16,8 +16,12 @@ #ifdef __NVCC__ #pragma push +#if (__CUDACC_VER_MAJOR__ >= 11) && (__CUDACC_VER_MINOR__ >= 5) +#pragma nv_diag_suppress declared_but_not_referenced // suppress "function was declared but never referenced warning" +#else #pragma diag_suppress declared_but_not_referenced // suppress "function was declared but never referenced warning" #endif +#endif #include "pugixml.h" diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index c7d68d73..3e06aa26 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -68,9 +68,16 @@ public: /////////////////////////////////////////////////////////////// // Support for MADWF tricks /////////////////////////////////////////////////////////////// - RealD Mass(void) { return mass; }; + RealD Mass(void) { return (mass_plus + mass_minus) / 2.0; }; + RealD MassPlus(void) { return mass_plus; }; + RealD MassMinus(void) { return mass_minus; }; void SetMass(RealD _mass) { - mass=_mass; + mass_plus=mass_minus=_mass; + SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs + } ; + void SetMass(RealD _mass_plus, RealD _mass_minus) { + mass_plus=_mass_plus; + mass_minus=_mass_minus; SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c); // Reset coeffs } ; void P(const FermionField &psi, FermionField &chi); @@ -108,7 +115,7 @@ public: void MeooeDag5D (const FermionField &in, FermionField &out); // protected: - RealD mass; + RealD mass_plus, mass_minus; // Save arguments to SetCoefficientsInternal Vector _gamma; diff --git a/Grid/qcd/action/fermion/CloverHelpers.h b/Grid/qcd/action/fermion/CloverHelpers.h new file mode 100644 index 00000000..57e71998 --- /dev/null +++ b/Grid/qcd/action/fermion/CloverHelpers.h @@ -0,0 +1,435 @@ +/************************************************************************************* + + 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); + + 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); + } + + EigenInvCloverOp = EigenCloverOp.inverse(); + 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); + pokeLocalSite(Qxinv, CTIv, lcoor); + }); + } + } + + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) { + return Helpers::Cmunu(U, lambda, mu, nu); + } + +}; + + +////////////////////////////////// +// 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; + + // Can this be avoided? + static void IdentityTimesC(const CloverField& in, RealD c) { + 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); + + Clover *= (1.0/diag_mass); + + // Taylor expansion, slow but generic + // Horner scheme: a0 + a1 x + a2 x^2 + .. = a0 + x (a1 + x(...)) + // qN = cN + // qn = cn + qn+1 X + std::vector cn(NMAX+1); + cn[0] = 1.0; + for (int i=1; i<=NMAX; i++) + cn[i] = cn[i-1] / RealD(i); + + ExpClover = Zero(); + IdentityTimesC(ExpClover, cn[NMAX]); + for (int i=NMAX-1; i>=0; i--) + ExpClover = ExpClover * Clover + cn[i]; + + // prepare inverse + CloverInv = (-1.0)*Clover; + + Clover = ExpClover * diag_mass; + + ExpClover = Zero(); + IdentityTimesC(ExpClover, cn[NMAX]); + for (int i=NMAX-1; i>=0; i--) + ExpClover = ExpClover * CloverInv + cn[i]; + + CloverInv = ExpClover * (1.0/diag_mass); + + } + + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) { + assert(0); + return lambda; + } + +}; + + +////////////////////////////////// +// 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 Exponentiate_Clover(CloverDiagonalField& Diagonal, + CloverTriangleField& Triangle, + RealD csw_t, RealD diag_mass) { + + // Do nothing + } + + // 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 MassTerm(CloverField& Clover, RealD diag_mass) { + // do nothing! + // mass term is multiplied to exp(Clover) below + } + + static int getNMAX(RealD prec, RealD R) { + /* compute stop condition for exponential */ + int NMAX=1; + RealD cond=R*R/2.; + + while (cond*std::exp(R)>prec) { + 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 ExponentiateHermitean6by6(const iMatrix &arg, const RealD& alpha, const std::vector& cN, const int Niter, iMatrix& dest){ + + typedef iMatrix mat; + + RealD qn[6]; + RealD qnold[6]; + RealD p[5]; + RealD trA2, trA3, trA4; + + mat A2, A3, A4, A5; + A2 = alpha * alpha * arg * arg; + A3 = alpha * arg * A2; + A4 = A2 * A2; + A5 = A2 * A3; + + trA2 = toReal( trace(A2) ); + trA3 = toReal( trace(A3) ); + trA4 = toReal( trace(A4)); + + p[0] = toReal( trace(A3 * A3)) / 6.0 - 0.125 * trA4 * trA2 - trA3 * trA3 / 18.0 + trA2 * trA2 * trA2/ 48.0; + p[1] = toReal( trace(A5)) / 5.0 - trA3 * trA2 / 6.0; + p[2] = toReal( trace(A4)) / 4.0 - 0.125 * trA2 * trA2; + p[3] = trA3 / 3.0; + p[4] = 0.5 * trA2; + + qnold[0] = cN[Niter]; + qnold[1] = 0.0; + qnold[2] = 0.0; + qnold[3] = 0.0; + qnold[4] = 0.0; + qnold[5] = 0.0; + + for(int i = Niter-1; i >= 0; i--) + { + qn[0] = p[0] * qnold[5] + cN[i]; + qn[1] = p[1] * qnold[5] + qnold[0]; + qn[2] = p[2] * qnold[5] + qnold[1]; + qn[3] = p[3] * qnold[5] + qnold[2]; + qn[4] = p[4] * qnold[5] + qnold[3]; + qn[5] = qnold[4]; + + qnold[0] = qn[0]; + qnold[1] = qn[1]; + qnold[2] = qn[2]; + qnold[3] = qn[3]; + qnold[4] = qn[4]; + qnold[5] = qn[5]; + } + + mat unit(1.0); + + dest = (qn[0] * unit + qn[1] * alpha * arg + qn[2] * A2 + qn[3] * A3 + qn[4] * A4 + qn[5] * A5); + + } + + static void Exponentiate_Clover(CloverDiagonalField& Diagonal, CloverTriangleField& Triangle, RealD csw_t, RealD diag_mass) { + + GridBase* grid = Diagonal.Grid(); + int NMAX = getNMAX(Diagonal, 3.*csw_t/diag_mass); + + // + // Implementation completely in Daniel's layout + // + + // Taylor expansion with Cayley-Hamilton recursion + // underlying Horner scheme as above + std::vector cn(NMAX+1); + cn[0] = 1.0; + for (int i=1; i<=NMAX; i++){ + cn[i] = cn[i-1] / RealD(i); + } + + // Taken over from Daniel's implementation + conformable(Diagonal, Triangle); + + long lsites = grid->lSites(); + { + typedef typename SiteCloverDiagonal::scalar_object scalar_object_diagonal; + typedef typename SiteCloverTriangle::scalar_object scalar_object_triangle; + typedef iMatrix mat; + + autoView(diagonal_v, Diagonal, CpuRead); + autoView(triangle_v, Triangle, CpuRead); + autoView(diagonalExp_v, Diagonal, CpuWrite); + autoView(triangleExp_v, Triangle, CpuWrite); + + thread_for(site, lsites, { // NOTE: Not on GPU because of (peek/poke)LocalSite + + mat srcCloverOpUL(0.0); // upper left block + mat srcCloverOpLR(0.0); // lower right block + mat ExpCloverOp; + + scalar_object_diagonal diagonal_tmp = Zero(); + scalar_object_diagonal diagonal_exp_tmp = Zero(); + scalar_object_triangle triangle_tmp = Zero(); + scalar_object_triangle triangle_exp_tmp = Zero(); + + Coordinate lcoor; + grid->LocalIndexToLocalCoor(site, lcoor); + + peekLocalSite(diagonal_tmp, diagonal_v, lcoor); + peekLocalSite(triangle_tmp, triangle_v, lcoor); + + int block; + block = 0; + for(int i = 0; i < 6; i++){ + for(int j = 0; j < 6; j++){ + if (i == j){ + srcCloverOpUL(i,j) = static_cast(TensorRemove(diagonal_tmp()(block)(i))); + } + else{ + srcCloverOpUL(i,j) = static_cast(TensorRemove(CompactHelpers::triangle_elem(triangle_tmp, block, i, j))); + } + } + } + block = 1; + for(int i = 0; i < 6; i++){ + for(int j = 0; j < 6; j++){ + if (i == j){ + srcCloverOpLR(i,j) = static_cast(TensorRemove(diagonal_tmp()(block)(i))); + } + else{ + srcCloverOpLR(i,j) = static_cast(TensorRemove(CompactHelpers::triangle_elem(triangle_tmp, block, i, j))); + } + } + } + + // exp(Clover) + + ExponentiateHermitean6by6(srcCloverOpUL,1.0/diag_mass,cn,NMAX,ExpCloverOp); + + block = 0; + for(int i = 0; i < 6; i++){ + for(int j = 0; j < 6; j++){ + if (i == j){ + diagonal_exp_tmp()(block)(i) = ExpCloverOp(i,j); + } + else if(i < j){ + triangle_exp_tmp()(block)(CompactHelpers::triangle_index(i, j)) = ExpCloverOp(i,j); + } + } + } + + ExponentiateHermitean6by6(srcCloverOpLR,1.0/diag_mass,cn,NMAX,ExpCloverOp); + + block = 1; + for(int i = 0; i < 6; i++){ + for(int j = 0; j < 6; j++){ + if (i == j){ + diagonal_exp_tmp()(block)(i) = ExpCloverOp(i,j); + } + else if(i < j){ + triangle_exp_tmp()(block)(CompactHelpers::triangle_index(i, j)) = ExpCloverOp(i,j); + } + } + } + + pokeLocalSite(diagonal_exp_tmp, diagonalExp_v, lcoor); + pokeLocalSite(triangle_exp_tmp, triangleExp_v, lcoor); + }); + } + + Diagonal *= diag_mass; + Triangle *= diag_mass; + } + + + static GaugeLinkField Cmunu(std::vector &U, GaugeLinkField &lambda, int mu, int nu) { + assert(0); + return lambda; + } + +}; + + +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..223ff9dd 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -138,38 +138,52 @@ typedef WilsonTMFermion WilsonTMFermionF; typedef WilsonTMFermion WilsonTMFermionD; // Clover fermions -typedef WilsonCloverFermion WilsonCloverFermionR; -typedef WilsonCloverFermion WilsonCloverFermionF; -typedef WilsonCloverFermion WilsonCloverFermionD; +template using WilsonClover = WilsonCloverFermion>; +template using WilsonExpClover = WilsonCloverFermion>; -typedef WilsonCloverFermion WilsonCloverAdjFermionR; -typedef WilsonCloverFermion WilsonCloverAdjFermionF; -typedef WilsonCloverFermion WilsonCloverAdjFermionD; +typedef WilsonClover WilsonCloverFermionR; +typedef WilsonClover WilsonCloverFermionF; +typedef WilsonClover WilsonCloverFermionD; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionR; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionF; -typedef WilsonCloverFermion WilsonCloverTwoIndexSymmetricFermionD; +typedef WilsonExpClover WilsonExpCloverFermionR; +typedef WilsonExpClover WilsonExpCloverFermionF; +typedef WilsonExpClover WilsonExpCloverFermionD; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionR; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionF; -typedef WilsonCloverFermion WilsonCloverTwoIndexAntiSymmetricFermionD; +typedef WilsonClover WilsonCloverAdjFermionR; +typedef WilsonClover WilsonCloverAdjFermionF; +typedef WilsonClover WilsonCloverAdjFermionD; + +typedef WilsonClover WilsonCloverTwoIndexSymmetricFermionR; +typedef WilsonClover WilsonCloverTwoIndexSymmetricFermionF; +typedef WilsonClover WilsonCloverTwoIndexSymmetricFermionD; + +typedef WilsonClover WilsonCloverTwoIndexAntiSymmetricFermionR; +typedef WilsonClover WilsonCloverTwoIndexAntiSymmetricFermionF; +typedef WilsonClover WilsonCloverTwoIndexAntiSymmetricFermionD; // Compact Clover fermions -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverFermionD; +template using CompactWilsonClover = CompactWilsonCloverFermion>; +template using CompactWilsonExpClover = CompactWilsonCloverFermion>; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverAdjFermionD; +typedef CompactWilsonClover CompactWilsonCloverFermionR; +typedef CompactWilsonClover CompactWilsonCloverFermionF; +typedef CompactWilsonClover CompactWilsonCloverFermionD; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexSymmetricFermionD; +typedef CompactWilsonExpClover CompactWilsonExpCloverFermionR; +typedef CompactWilsonExpClover CompactWilsonExpCloverFermionF; +typedef CompactWilsonExpClover CompactWilsonExpCloverFermionD; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionR; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionF; -typedef CompactWilsonCloverFermion CompactWilsonCloverTwoIndexAntiSymmetricFermionD; +typedef CompactWilsonClover CompactWilsonCloverAdjFermionR; +typedef CompactWilsonClover CompactWilsonCloverAdjFermionF; +typedef CompactWilsonClover CompactWilsonCloverAdjFermionD; + +typedef CompactWilsonClover CompactWilsonCloverTwoIndexSymmetricFermionR; +typedef CompactWilsonClover CompactWilsonCloverTwoIndexSymmetricFermionF; +typedef CompactWilsonClover CompactWilsonCloverTwoIndexSymmetricFermionD; + +typedef CompactWilsonClover CompactWilsonCloverTwoIndexAntiSymmetricFermionR; +typedef CompactWilsonClover CompactWilsonCloverTwoIndexAntiSymmetricFermionF; +typedef CompactWilsonClover 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/WilsonCloverTypes.h b/Grid/qcd/action/fermion/WilsonCloverTypes.h index a5a95d93..5ce3b91b 100644 --- a/Grid/qcd/action/fermion/WilsonCloverTypes.h +++ b/Grid/qcd/action/fermion/WilsonCloverTypes.h @@ -47,8 +47,6 @@ class CompactWilsonCloverTypes { public: INHERIT_IMPL_TYPES(Impl); - static_assert(Nd == 4 && Nc == 3 && Ns == 4 && Impl::Dimension == 3, "Wrong dimensions"); - static constexpr int Nred = Nc * Nhs; // 6 static constexpr int Nblock = Nhs; // 2 static constexpr int Ndiagonal = Nred; // 6 diff --git a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h index 1e2b13e3..51a7990c 100644 --- a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h @@ -47,7 +47,7 @@ CayleyFermion5D::CayleyFermion5D(GaugeField &_Umu, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,_M5,p), - mass(_mass) + mass_plus(_mass), mass_minus(_mass) { } @@ -209,8 +209,8 @@ void CayleyFermion5D::M5D (const FermionField &psi, FermionField &chi) { int Ls=this->Ls; Vector diag (Ls,1.0); - Vector upper(Ls,-1.0); upper[Ls-1]=mass; - Vector lower(Ls,-1.0); lower[0] =mass; + Vector upper(Ls,-1.0); upper[Ls-1]=mass_minus; + Vector lower(Ls,-1.0); lower[0] =mass_plus; M5D(psi,chi,chi,lower,diag,upper); } template @@ -220,8 +220,8 @@ void CayleyFermion5D::Meooe5D (const FermionField &psi, FermionField &D Vector diag = bs; Vector upper= cs; Vector lower= cs; - upper[Ls-1]=-mass*upper[Ls-1]; - lower[0] =-mass*lower[0]; + upper[Ls-1]=-mass_minus*upper[Ls-1]; + lower[0] =-mass_plus*lower[0]; M5D(psi,psi,Din,lower,diag,upper); } // FIXME Redunant with the above routine; check this and eliminate @@ -235,8 +235,8 @@ template void CayleyFermion5D::Meo5D (const FermionField & upper[i]=-ceo[i]; lower[i]=-ceo[i]; } - upper[Ls-1]=-mass*upper[Ls-1]; - lower[0] =-mass*lower[0]; + upper[Ls-1]=-mass_minus*upper[Ls-1]; + lower[0] =-mass_plus*lower[0]; M5D(psi,psi,chi,lower,diag,upper); } template @@ -250,8 +250,8 @@ void CayleyFermion5D::Mooee (const FermionField &psi, FermionField & upper[i]=-cee[i]; lower[i]=-cee[i]; } - upper[Ls-1]=-mass*upper[Ls-1]; - lower[0] =-mass*lower[0]; + upper[Ls-1]=-mass_minus*upper[Ls-1]; + lower[0] =-mass_plus*lower[0]; M5D(psi,psi,chi,lower,diag,upper); } template @@ -266,9 +266,9 @@ void CayleyFermion5D::MooeeDag (const FermionField &psi, FermionField & // Assemble the 5d matrix if ( s==0 ) { upper[s] = -cee[s+1] ; - lower[s] = mass*cee[Ls-1]; + lower[s] = mass_minus*cee[Ls-1]; } else if ( s==(Ls-1)) { - upper[s] = mass*cee[0]; + upper[s] = mass_plus*cee[0]; lower[s] = -cee[s-1]; } else { upper[s]=-cee[s+1]; @@ -291,8 +291,8 @@ void CayleyFermion5D::M5Ddag (const FermionField &psi, FermionField &chi) Vector diag(Ls,1.0); Vector upper(Ls,-1.0); Vector lower(Ls,-1.0); - upper[Ls-1]=-mass*upper[Ls-1]; - lower[0] =-mass*lower[0]; + upper[Ls-1]=-mass_plus*upper[Ls-1]; + lower[0] =-mass_minus*lower[0]; M5Ddag(psi,chi,chi,lower,diag,upper); } @@ -307,9 +307,9 @@ void CayleyFermion5D::MeooeDag5D (const FermionField &psi, FermionField for (int s=0;s::SetCoefficientsInternal(RealD zolo_hi,Vector::SetCoefficientsInternal(RealD zolo_hi,Vector::SetCoefficientsInternal(RealD zolo_hi,Vector::ContractConservedCurrent( PropagatorField &q_in_1, Current curr_type, unsigned int mu) { + + assert(mass_plus == mass_minus); + RealD mass = mass_plus; + #if (!defined(GRID_HIP)) Gamma::Algebra Gmu [] = { Gamma::Algebra::GammaX, @@ -777,6 +781,8 @@ void CayleyFermion5D::SeqConservedCurrent(PropagatorField &q_in, assert(mu>=0); assert(mu #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) @@ -58,50 +59,55 @@ CompactWilsonCloverFermion::CompactWilsonCloverFermion(GaugeField& _Umu, , BoundaryMask(&Fgrid) , BoundaryMaskEven(&Hgrid), BoundaryMaskOdd(&Hgrid) { + assert(Nd == 4 && Nc == 3 && Ns == 4 && Impl::Dimension == 3); + csw_r *= 0.5; csw_t *= 0.5; if (clover_anisotropy.isAnisotropic) csw_r /= clover_anisotropy.xi_0; ImportGauge(_Umu); - if (open_boundaries) + if (open_boundaries) { + this->BoundaryMaskEven.Checkerboard() = Even; + this->BoundaryMaskOdd.Checkerboard() = Odd; 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 +115,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 +124,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 +150,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 +169,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 +257,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 +267,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 +291,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,22 +324,27 @@ 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); - // Possible modify the boundary values + // Exponentiate the clover (nothing happens in case of the standard clover) double t5 = usecond(); + CloverHelpers::Exponentiate_Clover(Diagonal, Triangle, csw_t, this->diag_mass); + + // Possible modify the boundary values + double t6 = usecond(); if(open_boundaries) CompactHelpers::ModifyBoundaries(Diagonal, Triangle, csw_t, cF, this->diag_mass); - // Invert the clover term in the improved layout - double t6 = usecond(); + // Invert the Clover term (explicit inversion needed for the improvement in case of open boundary conditions) + double t7 = usecond(); CompactHelpers::Invert(Diagonal, Triangle, DiagonalInv, TriangleInv); // Fill the remaining clover fields - double t7 = usecond(); + double t8 = usecond(); pickCheckerboard(Even, DiagonalEven, Diagonal); pickCheckerboard(Even, TriangleEven, Triangle); pickCheckerboard(Odd, DiagonalOdd, Diagonal); @@ -344,20 +355,19 @@ void CompactWilsonCloverFermion::ImportGauge(const GaugeField& _Umu) { pickCheckerboard(Odd, TriangleInvOdd, TriangleInv); // Report timings - double t8 = usecond(); -#if 0 - std::cout << GridLogMessage << "CompactWilsonCloverFermion::ImportGauge timings:" - << " WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 - << ", allocations = " << (t2 - t1) / 1e6 - << ", field strength = " << (t3 - t2) / 1e6 - << ", fill clover = " << (t4 - t3) / 1e6 - << ", convert = " << (t5 - t4) / 1e6 - << ", boundaries = " << (t6 - t5) / 1e6 - << ", inversions = " << (t7 - t6) / 1e6 - << ", pick cbs = " << (t8 - t7) / 1e6 - << ", total = " << (t8 - t0) / 1e6 - << std::endl; -#endif + double t9 = usecond(); + + std::cout << GridLogDebug << "CompactWilsonCloverFermion::ImportGauge timings:" << std::endl; + std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl; + std::cout << GridLogDebug << "allocations = " << (t2 - t1) / 1e6 << std::endl; + std::cout << GridLogDebug << "field strength = " << (t3 - t2) / 1e6 << std::endl; + std::cout << GridLogDebug << "fill clover = " << (t4 - t3) / 1e6 << std::endl; + std::cout << GridLogDebug << "convert = " << (t5 - t4) / 1e6 << std::endl; + std::cout << GridLogDebug << "exponentiation = " << (t6 - t5) / 1e6 << std::endl; + std::cout << GridLogDebug << "boundaries = " << (t7 - t6) / 1e6 << std::endl; + std::cout << GridLogDebug << "inversions = " << (t8 - t7) / 1e6 << std::endl; + std::cout << GridLogDebug << "pick cbs = " << (t9 - t8) / 1e6 << std::endl; + std::cout << GridLogDebug << "total = " << (t9 - t0) / 1e6 << std::endl; } NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h index e4d2e736..2a7e7535 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,11 @@ 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); 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,48 +148,44 @@ 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:" - << " 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 + std::cout << GridLogDebug << "WilsonCloverFermion::ImportGauge timings:" << std::endl; + std::cout << GridLogDebug << "WilsonFermion::Importgauge = " << (t1 - t0) / 1e6 << std::endl; + std::cout << GridLogDebug << "allocations = " << (t2 - t1) / 1e6 << std::endl; + std::cout << GridLogDebug << "field strength = " << (t3 - t2) / 1e6 << std::endl; + std::cout << GridLogDebug << "fill clover = " << (t4 - t3) / 1e6 << std::endl; + std::cout << GridLogDebug << "instantiation = " << (t5 - t4) / 1e6 << std::endl; + std::cout << GridLogDebug << "pick cbs = " << (t6 - t5) / 1e6 << std::endl; + std::cout << GridLogDebug << "total = " << (t6 - t0) / 1e6 << std::endl; } -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 +238,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 +309,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 +320,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/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h index 84ac25c1..c958019d 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h @@ -4,12 +4,13 @@ Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/qcd/action/fermion/WilsonFermion.cc -Copyright (C) 2015 +Copyright (C) 2022 Author: Peter Boyle Author: Peter Boyle Author: Peter Boyle Author: paboyle +Author: Fabian Joswig 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 @@ -599,11 +600,47 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, Current curr_type, unsigned int mu) { + if(curr_type != Current::Vector) + { + std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl; + exit(1); + } + Gamma g5(Gamma::Algebra::Gamma5); conformable(_grid, q_in_1.Grid()); conformable(_grid, q_in_2.Grid()); conformable(_grid, q_out.Grid()); - assert(0); + auto UGrid= this->GaugeGrid(); + + PropagatorField tmp_shifted(UGrid); + PropagatorField g5Lg5(UGrid); + PropagatorField R(UGrid); + PropagatorField gmuR(UGrid); + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT, + }; + Gamma gmu=Gamma(Gmu[mu]); + + g5Lg5=g5*q_in_1*g5; + tmp_shifted=Cshift(q_in_2,mu,1); + Impl::multLinkField(R,this->Umu,tmp_shifted,mu); + gmuR=gmu*R; + + q_out=adj(g5Lg5)*R; + q_out-=adj(g5Lg5)*gmuR; + + tmp_shifted=Cshift(q_in_1,mu,1); + Impl::multLinkField(g5Lg5,this->Umu,tmp_shifted,mu); + g5Lg5=g5*g5Lg5*g5; + R=q_in_2; + gmuR=gmu*R; + + q_out-=adj(g5Lg5)*R; + q_out-=adj(g5Lg5)*gmuR; } @@ -617,9 +654,51 @@ void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, unsigned int tmax, ComplexField &lattice_cmplx) { + if(curr_type != Current::Vector) + { + std::cout << GridLogError << "Only the conserved vector current is implemented so far." << std::endl; + exit(1); + } + + int tshift = (mu == Nd-1) ? 1 : 0; + unsigned int LLt = GridDefaultLatt()[Tp]; conformable(_grid, q_in.Grid()); conformable(_grid, q_out.Grid()); - assert(0); + auto UGrid= this->GaugeGrid(); + + PropagatorField tmp(UGrid); + PropagatorField Utmp(UGrid); + PropagatorField L(UGrid); + PropagatorField zz (UGrid); + zz=Zero(); + LatticeInteger lcoor(UGrid); LatticeCoordinate(lcoor,Nd-1); + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT, + }; + Gamma gmu=Gamma(Gmu[mu]); + + tmp = Cshift(q_in,mu,1); + Impl::multLinkField(Utmp,this->Umu,tmp,mu); + tmp = ( Utmp*lattice_cmplx - gmu*Utmp*lattice_cmplx ); // Forward hop + tmp = where((lcoor>=tmin),tmp,zz); // Mask the time + q_out = where((lcoor<=tmax),tmp,zz); // Position of current complicated + + tmp = q_in *lattice_cmplx; + tmp = Cshift(tmp,mu,-1); + Impl::multLinkField(Utmp,this->Umu,tmp,mu+Nd); // Adjoint link + tmp = -( Utmp + gmu*Utmp ); + // Mask the time + if (tmax == LLt - 1 && tshift == 1){ // quick fix to include timeslice 0 if tmax + tshift is over the last timeslice + unsigned int t0 = 0; + tmp = where(((lcoor==t0) || (lcoor>=tmin+tshift)),tmp,zz); + } else { + tmp = where((lcoor>=tmin+tshift),tmp,zz); + } + q_out+= where((lcoor<=tmax+tshift),tmp,zz); // Position of current complicated } NAMESPACE_END(Grid); 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/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplD/WilsonKernelsInstantiationWilsonAdjImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonAdjImplF/WilsonKernelsInstantiationWilsonAdjImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file 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); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonImplD/WilsonKernelsInstantiationWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonImplF/WilsonKernelsInstantiationWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexAntiSymmetricImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/WilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationWilsonTwoIndexSymmetricImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc b/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/ZWilsonImplD/WilsonKernelsInstantiationZWilsonImplD.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc deleted file mode 100644 index f0b15e3b..00000000 --- a/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc +++ /dev/null @@ -1,51 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: ./lib/qcd/action/fermion/WilsonKernels.cc - -Copyright (C) 2015, 2020 - -Author: Peter Boyle -Author: Peter Boyle -Author: paboyle -Author: Nils Meyer Regensburg University - -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 - -#ifndef AVX512 -#ifndef QPX -#ifndef A64FX -#ifndef A64FXFIXEDSIZE -#include -#endif -#endif -#endif -#endif - -NAMESPACE_BEGIN(Grid); - -#include "impl.h" -template class WilsonKernels; - -NAMESPACE_END(Grid); diff --git a/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc b/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc new file mode 120000 index 00000000..01c35e7b --- /dev/null +++ b/Grid/qcd/action/fermion/instantiation/ZWilsonImplF/WilsonKernelsInstantiationZWilsonImplF.cc @@ -0,0 +1 @@ +../WilsonKernelsInstantiation.cc.master \ No newline at end of file diff --git a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh index a09fd420..d6845e75 100755 --- a/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh +++ b/Grid/qcd/action/fermion/instantiation/generate_instantiations.sh @@ -18,6 +18,10 @@ WILSON_IMPL_LIST=" \ GparityWilsonImplF \ GparityWilsonImplD " +COMPACT_WILSON_IMPL_LIST=" \ + WilsonImplF \ + WilsonImplD " + DWF_IMPL_LIST=" \ WilsonImplF \ WilsonImplD \ @@ -40,13 +44,23 @@ EOF done -CC_LIST="WilsonCloverFermionInstantiation CompactWilsonCloverFermionInstantiation WilsonFermionInstantiation WilsonKernelsInstantiation WilsonTMFermionInstantiation" +CC_LIST="WilsonCloverFermionInstantiation WilsonFermionInstantiation WilsonKernelsInstantiation WilsonTMFermionInstantiation" for impl in $WILSON_IMPL_LIST do for f in $CC_LIST do - ln -f -s ../$f.cc.master $impl/$f$impl.cc + ln -f -s ../$f.cc.master $impl/$f$impl.cc +done +done + +CC_LIST="CompactWilsonCloverFermionInstantiation" + +for impl in $COMPACT_WILSON_IMPL_LIST +do +for f in $CC_LIST +do + ln -f -s ../$f.cc.master $impl/$f$impl.cc done done @@ -63,14 +77,14 @@ for impl in $DWF_IMPL_LIST $GDWF_IMPL_LIST do for f in $CC_LIST do - ln -f -s ../$f.cc.master $impl/$f$impl.cc + ln -f -s ../$f.cc.master $impl/$f$impl.cc done done # overwrite the .cc file in Gparity directories for impl in $GDWF_IMPL_LIST do - ln -f -s ../WilsonKernelsInstantiationGparity.cc.master $impl/WilsonKernelsInstantiation$impl.cc + ln -f -s ../WilsonKernelsInstantiationGparity.cc.master $impl/WilsonKernelsInstantiation$impl.cc done @@ -84,7 +98,7 @@ for impl in $STAG_IMPL_LIST do for f in $CC_LIST do - ln -f -s ../$f.cc.master $impl/$f$impl.cc + ln -f -s ../$f.cc.master $impl/$f$impl.cc done done diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index 2b3384da..29184a88 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -55,12 +55,12 @@ public: } } - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { GridBase *grid = Umu.Grid(); GaugeMat xform(grid); - SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog); + SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier,orthog,err_on_no_converge); } - static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1) { + static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false,int orthog=-1,bool err_on_no_converge=true) { GridBase *grid = Umu.Grid(); @@ -122,6 +122,8 @@ public: } } + std::cout << GridLogError << "Gauge fixing did not converge in " << maxiter << " iterations." << std::endl; + if (err_on_no_converge) assert(0); }; static Real SteepestDescentStep(std::vector &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) { GridBase *grid = U[0].Grid(); diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 0367c9fa..e002e3d5 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -125,7 +125,6 @@ public: return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } - ////////////////////////////////////////////////// // average over all x,y,z the temporal loop ////////////////////////////////////////////////// @@ -165,7 +164,7 @@ public: double vol = Umu.Grid()->gSites(); - return p.real() / vol / 4.0 / 3.0; + return p.real() / vol / (4.0 * Nc ) ; }; ////////////////////////////////////////////////// diff --git a/benchmarks/Benchmark_mooee.cc b/benchmarks/Benchmark_mooee.cc index 0aaccecc..289af41d 100644 --- a/benchmarks/Benchmark_mooee.cc +++ b/benchmarks/Benchmark_mooee.cc @@ -81,8 +81,8 @@ int main (int argc, char ** argv) Vector diag = Dw.bs; Vector upper= Dw.cs; Vector lower= Dw.cs; - upper[Ls-1]=-Dw.mass*upper[Ls-1]; - lower[0] =-Dw.mass*lower[0]; + upper[Ls-1]=-Dw.mass_minus*upper[Ls-1]; + lower[0] =-Dw.mass_plus*lower[0]; LatticeFermion r_eo(FGrid); LatticeFermion src_e (FrbGrid); diff --git a/benchmarks/Benchmark_wilson_sweep.cc b/benchmarks/Benchmark_wilson_sweep.cc index 986a1831..50ad1d03 100644 --- a/benchmarks/Benchmark_wilson_sweep.cc +++ b/benchmarks/Benchmark_wilson_sweep.cc @@ -44,6 +44,13 @@ void bench_wilson ( double const volume, int const dag ); +void bench_wilson_eo ( + LatticeFermion & src, + LatticeFermion & result, + WilsonFermionR & Dw, + double const volume, + int const dag ); + int main (int argc, char ** argv) { Grid_init(&argc,&argv); @@ -110,8 +117,8 @@ int main (int argc, char ** argv) bench_wilson(src,result,Dw,volume,DaggerYes); std::cout << "\t"; // EO - bench_wilson(src,result,Dw,volume,DaggerNo); - bench_wilson(src,result,Dw,volume,DaggerYes); + bench_wilson_eo(src_o,result_e,Dw,volume,DaggerNo); + bench_wilson_eo(src_o,result_e,Dw,volume,DaggerYes); std::cout << std::endl; } } diff --git a/configure.ac b/configure.ac index 406b0b74..9ab0595a 100644 --- a/configure.ac +++ b/configure.ac @@ -159,7 +159,7 @@ case ${ac_ZMOBIUS} in esac ############### Nc AC_ARG_ENABLE([Nc], - [AC_HELP_STRING([--enable-Nc=2|3|4], [enable number of colours])], + [AC_HELP_STRING([--enable-Nc=2|3|4|5], [enable number of colours])], [ac_Nc=${enable_Nc}], [ac_Nc=3]) case ${ac_Nc} in diff --git a/tests/IO/Test_nersc_io.cc b/tests/IO/Test_nersc_io.cc index c15c320e..ae5b9c0d 100644 --- a/tests/IO/Test_nersc_io.cc +++ b/tests/IO/Test_nersc_io.cc @@ -147,7 +147,7 @@ int main (int argc, char ** argv) Complex p = TensorRemove(Tp); std::cout< WImpl; - typedef WilsonCloverFermion WilsonCloverOperator; - typedef CompactWilsonCloverFermion CompactWilsonCloverOperator; + typedef WilsonCloverFermion> WilsonCloverOperator; + typedef CompactWilsonCloverFermion> CompactWilsonCloverOperator; typedef typename WilsonCloverOperator::FermionField Fermion; typedef typename WilsonCloverOperator::GaugeField Gauge; diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index e044378c..8f18ed04 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -9,6 +9,7 @@ Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle Author: Guido Cossu +Author: Jamie Hudspith 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 @@ -42,14 +43,14 @@ directory using namespace std; using namespace Grid; - ; +; int main(int argc, char** argv) { Grid_init(&argc, &argv); std::vector latt({4, 4, 4, 8}); GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( - latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); + latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); @@ -60,15 +61,19 @@ int main(int argc, char** argv) { << std::endl; SU2::printGenerators(); std::cout << "Dimension of adjoint representation: "<< SU2Adjoint::Dimension << std::endl; + + // guard as this code fails to compile for Nc != 3 +#if (Nc == 3) + SU2Adjoint::printGenerators(); SU2::testGenerators(); SU2Adjoint::testGenerators(); - + std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; SU3::printGenerators(); std::cout << "Dimension of adjoint representation: "<< SU3Adjoint::Dimension << std::endl; SU3Adjoint::printGenerators(); @@ -111,12 +116,10 @@ int main(int argc, char** argv) { // AdjointRepresentation has the predefined number of colours Nc // Representations RepresentationTypes(grid); - - LatticeGaugeField U(grid), V(grid); SU3::HotConfiguration(gridRNG, U); SU3::HotConfiguration(gridRNG, V); - + // Adjoint representation // Test group structure // (U_f * V_f)_r = U_r * V_r @@ -127,17 +130,17 @@ int main(int argc, char** argv) { SU3::LatticeMatrix Vmu = peekLorentz(V,mu); pokeLorentz(UV,Umu*Vmu, mu); } - + AdjRep.update_representation(UV); typename AdjointRep::LatticeField UVr = AdjRep.U; // (U_f * V_f)_r - - + + AdjRep.update_representation(U); typename AdjointRep::LatticeField Ur = AdjRep.U; // U_r - + AdjRep.update_representation(V); typename AdjointRep::LatticeField Vr = AdjRep.U; // V_r - + typename AdjointRep::LatticeField UrVr(grid); UrVr = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -145,10 +148,10 @@ int main(int argc, char** argv) { typename AdjointRep::LatticeMatrix Vrmu = peekLorentz(Vr,mu); pokeLorentz(UrVr,Urmu*Vrmu, mu); } - + typename AdjointRep::LatticeField Diff_check = UVr - UrVr; std::cout << GridLogMessage << "Group structure SU("<::AdjointLieAlgebraMatrix(h_adj,Ar); - + // Re-extract h_adj SU3::LatticeAlgebraVector h_adj2(grid); SU_Adjoint::projectOnAlgebra(h_adj2, Ar); SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2; std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl; - + // Exponentiate typename AdjointRep::LatticeMatrix Uadj(grid); Uadj = expMat(Ar, 1.0, 16); - + typename AdjointRep::LatticeMatrix uno(grid); uno = 1.0; // Check matrix Uadj, must be real orthogonal typename AdjointRep::LatticeMatrix Ucheck = Uadj - conjugate(Uadj); std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck) - << std::endl; - + << std::endl; + Ucheck = Uadj * adj(Uadj) - uno; std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck) - << std::endl; + << std::endl; Ucheck = adj(Uadj) * Uadj - uno; std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck) - << std::endl; - - + << std::endl; + // Construct the fundamental matrix in the group SU3::LatticeMatrix Af(grid); SU3::FundamentalLieAlgebraMatrix(h_adj,Af); @@ -193,72 +195,65 @@ int main(int argc, char** argv) { SU3::LatticeMatrix UnitCheck(grid); UnitCheck = Ufund * adj(Ufund) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) - << std::endl; + << std::endl; UnitCheck = adj(Ufund) * Ufund - uno_f; std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck) - << std::endl; - + << std::endl; // Tranform to the adjoint representation U = Zero(); // fill this with only one direction pokeLorentz(U,Ufund,0); // the representation transf acts on full gauge fields - + AdjRep.update_representation(U); Ur = AdjRep.U; // U_r typename AdjointRep::LatticeMatrix Ur0 = peekLorentz(Ur,0); // this should be the same as Uadj - + typename AdjointRep::LatticeMatrix Diff_check_mat = Ur0 - Uadj; std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl; - - - // TwoIndexRep tests - std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* eS^{ij} base for SU(2)" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Dimension of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::printBase(); - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; + std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::printGenerators(); - std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl; + std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl; SU2TwoIndexSymm::testGenerators(); std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "* eAS^{ij} base for SU(2)" << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::printBase(); - std::cout << GridLogMessage << "*********************************************" - << std::endl; - std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; + std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::printGenerators(); std::cout << GridLogMessage << "Test of Two Index anti-Symmetric Generators: "<< SU2TwoIndexAntiSymm::Dimension << std::endl; SU2TwoIndexAntiSymm::testGenerators(); - - std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors" - << std::endl; + << std::endl; // Projectors SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid); random(gridRNG,Gauss2); @@ -276,13 +271,13 @@ int main(int argc, char** argv) { SU3::LatticeAlgebraVector diff2 = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; - std::cout << GridLogMessage << "*********************************************" - << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors" - << std::endl; + << std::endl; // Projectors SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid); random(gridRNG,Gauss2a); @@ -300,11 +295,11 @@ int main(int argc, char** argv) { SU3::LatticeAlgebraVector diff2a = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl; std::cout << GridLogMessage << "*********************************************" - << std::endl; + << std::endl; std::cout << GridLogMessage << "Two index Symmetric: Checking Group Structure" - << std::endl; + << std::endl; // Testing HMC representation classes TwoIndexRep< Nc, Symmetric > TIndexRep(grid); @@ -313,7 +308,7 @@ int main(int argc, char** argv) { LatticeGaugeField U2(grid), V2(grid); SU3::HotConfiguration(gridRNG, U2); SU3::HotConfiguration(gridRNG, V2); - + LatticeGaugeField UV2(grid); UV2 = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -321,16 +316,16 @@ int main(int argc, char** argv) { SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu); pokeLorentz(UV2,Umu2*Vmu2, mu); } - + TIndexRep.update_representation(UV2); typename TwoIndexRep< Nc, Symmetric >::LatticeField UVr2 = TIndexRep.U; // (U_f * V_f)_r - + TIndexRep.update_representation(U2); typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2 = TIndexRep.U; // U_r - + TIndexRep.update_representation(V2); typename TwoIndexRep< Nc, Symmetric >::LatticeField Vr2 = TIndexRep.U; // V_r - + typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2Vr2(grid); Ur2Vr2 = Zero(); for (int mu = 0; mu < Nd; mu++) { @@ -338,11 +333,11 @@ int main(int argc, char** argv) { typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Vrmu2 = peekLorentz(Vr2,mu); pokeLorentz(Ur2Vr2,Urmu2*Vrmu2, mu); } - + typename TwoIndexRep< Nc, Symmetric >::LatticeField Diff_check2 = UVr2 - Ur2Vr2; std::cout << GridLogMessage << "Group structure SU("<::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym); - + // Re-extract h_sym SU3::LatticeAlgebraVector h_sym2(grid); SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym); SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl; - - + // Exponentiate typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix U2iS(grid); U2iS = expMat(Ar_sym, 1.0, 16); - + typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix uno2iS(grid); uno2iS = 1.0; // Check matrix U2iS, must be real orthogonal typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ucheck2iS = U2iS - conjugate(U2iS); std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iS) - << std::endl; - + << std::endl; + Ucheck2iS = U2iS * adj(U2iS) - uno2iS; std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iS) - << std::endl; + << std::endl; Ucheck2iS = adj(U2iS) * U2iS - uno2iS; std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iS) - << std::endl; - - - + << std::endl; + // Construct the fundamental matrix in the group SU3::LatticeMatrix Af_sym(grid); SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym); @@ -386,147 +378,137 @@ int main(int argc, char** argv) { SU3::LatticeMatrix UnitCheck2(grid); UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2) - << std::endl; + << std::endl; UnitCheck2 = adj(Ufund2) * Ufund2 - uno_f; std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2) - << std::endl; - - + << std::endl; + + // Tranform to the 2Index Sym representation U = Zero(); // fill this with only one direction pokeLorentz(U,Ufund2,0); // the representation transf acts on full gauge fields - + TIndexRep.update_representation(U); Ur2 = TIndexRep.U; // U_r typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ur02 = peekLorentz(Ur2,0); // this should be the same as U2iS - + typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Diff_check_mat2 = Ur02 - U2iS; std::cout << GridLogMessage << "Projections structure check group difference (Two Index Symmetric): " << norm2(Diff_check_mat2) << std::endl; - - - - + if (TwoIndexRep::Dimension != 1){ - std::cout << GridLogMessage << "*********************************************" - << std::endl; + std::cout << GridLogMessage << "*********************************************" + << std::endl; - std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure" - << std::endl; - // Testing HMC representation classes - TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid); + std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure" + << std::endl; + // Testing HMC representation classes + TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid); - // Test group structure - // (U_f * V_f)_r = U_r * V_r - LatticeGaugeField U2A(grid), V2A(grid); - SU3::HotConfiguration(gridRNG, U2A); - SU3::HotConfiguration(gridRNG, V2A); + // Test group structure + // (U_f * V_f)_r = U_r * V_r + LatticeGaugeField U2A(grid), V2A(grid); + SU3::HotConfiguration(gridRNG, U2A); + SU3::HotConfiguration(gridRNG, V2A); - LatticeGaugeField UV2A(grid); - UV2A = Zero(); - for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); - pokeLorentz(UV2A,Umu2A*Vmu2A, mu); + LatticeGaugeField UV2A(grid); + UV2A = Zero(); + for (int mu = 0; mu < Nd; mu++) { + SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); + SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); + pokeLorentz(UV2A,Umu2A*Vmu2A, mu); + } + + TIndexRep.update_representation(UV2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r + + TIndexRep.update_representation(U2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r + + TIndexRep.update_representation(V2A); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r + + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid); + Ur2Vr2A = Zero(); + for (int mu = 0; mu < Nd; mu++) { + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu); + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu); + pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu); + } + + typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A; + std::cout << GridLogMessage << "Group structure SU("<::LatticeMatrix Ar_Asym(grid); + random(gridRNG,h_Asym); + h_Asym = real(h_Asym); + SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); + + // Re-extract h_sym + SU3::LatticeAlgebraVector h_Asym2(grid); + SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); + SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; + std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; + + + // Exponentiate + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid); + U2iAS = expMat(Ar_Asym, 1.0, 16); + + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid); + uno2iAS = 1.0; + // Check matrix U2iS, must be real orthogonal + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS); + std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS) + << std::endl; + + Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS; + std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS) + << std::endl; + Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS; + std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS) + << std::endl; + + + + // Construct the fundamental matrix in the group + SU3::LatticeMatrix Af_Asym(grid); + SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); + SU3::LatticeMatrix Ufund2A(grid); + Ufund2A = expMat(Af_Asym, 1.0, 16); + SU3::LatticeMatrix UnitCheck2A(grid); + UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; + std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) + << std::endl; + UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f; + std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A) + << std::endl; + + + // Tranform to the 2Index Sym representation + U = Zero(); // fill this with only one direction + pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields + + TIndexRepA.update_representation(U); + Ur2A = TIndexRepA.U; // U_r + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS + + typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS; + std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl; + + } else { + std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests " + "because representation is trivial (dim = 1)" + << std::endl; } - - TIndexRep.update_representation(UV2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U; // (U_f * V_f)_r - - TIndexRep.update_representation(U2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U; // U_r - - TIndexRep.update_representation(V2A); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U; // V_r - - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid); - Ur2Vr2A = Zero(); - for (int mu = 0; mu < Nd; mu++) { - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu); - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu); - pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu); - } - - typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A; - std::cout << GridLogMessage << "Group structure SU("<::LatticeMatrix Ar_Asym(grid); - random(gridRNG,h_Asym); - h_Asym = real(h_Asym); - SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym); - - // Re-extract h_sym - SU3::LatticeAlgebraVector h_Asym2(grid); - SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); - SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; - std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl; - - - // Exponentiate - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid); - U2iAS = expMat(Ar_Asym, 1.0, 16); - - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid); - uno2iAS = 1.0; - // Check matrix U2iS, must be real orthogonal - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS); - std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS) - << std::endl; - - Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS; - std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS) - << std::endl; - Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS; - std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS) - << std::endl; - - - - // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_Asym(grid); - SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); - SU3::LatticeMatrix Ufund2A(grid); - Ufund2A = expMat(Af_Asym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2A(grid); - UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; - std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) - << std::endl; - UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f; - std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A) - << std::endl; - - - // Tranform to the 2Index Sym representation - U = Zero(); // fill this with only one direction - pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields - - TIndexRepA.update_representation(U); - Ur2A = TIndexRepA.U; // U_r - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS - - typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS; - std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl; - -} else { - std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests " - "because representation is trivial (dim = 1)" - << std::endl; -} - - - - - - - - +#endif Grid_finalize(); } diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index 6644be1a..10423159 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -122,14 +122,15 @@ int main (int argc, char ** argv) std::cout << "Determinant defect before projection " < + Fabian Joswig 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 @@ -67,8 +68,6 @@ int main(int argc, char **argv) tmp = Zero(); FermionField err(&Grid); err = Zero(); - FermionField err2(&Grid); - err2 = Zero(); FermionField phi(&Grid); random(pRNG, phi); FermionField chi(&Grid); @@ -77,6 +76,8 @@ int main(int argc, char **argv) SU::HotConfiguration(pRNG, Umu); std::vector U(4, &Grid); + double tolerance = 1e-4; + double volume = 1; for (int mu = 0; mu < Nd; mu++) { @@ -88,7 +89,7 @@ int main(int argc, char **argv) RealD csw_t = 1.0; WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params); - //Dwc.ImportGauge(Umu); // not necessary, included in the constructor + CompactWilsonCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params); std::cout << GridLogMessage << "==========================================================" << std::endl; std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl; @@ -112,7 +113,24 @@ int main(int argc, char **argv) setCheckerboard(r_eo, r_e); err = ref - r_eo; - std::cout << GridLogMessage << "EO norm diff " << norm2(err) << " " << norm2(ref) << " " << norm2(r_eo) << std::endl; + std::cout << GridLogMessage << "EO norm diff\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl; + assert(fabs(norm2(err)) < tolerance); + + + + Dwc_compact.Meooe(src_e, r_o); + std::cout << GridLogMessage << "Applied Meo" << std::endl; + Dwc_compact.Meooe(src_o, r_e); + std::cout << GridLogMessage << "Applied Moe" << std::endl; + Dwc_compact.Dhop(src, ref, DaggerNo); + + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff compact\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl; + assert(fabs(norm2(err)) < tolerance); + std::cout << GridLogMessage << "==============================================================" << std::endl; std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; @@ -152,6 +170,22 @@ int main(int argc, char **argv) std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl; std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl; + Dwc_compact.Meooe(chi_e, dchi_o); + Dwc_compact.Meooe(chi_o, dchi_e); + Dwc_compact.MeooeDag(phi_e, dphi_o); + Dwc_compact.MeooeDag(phi_o, dphi_e); + + pDce = innerProduct(phi_e, dchi_e); + pDco = innerProduct(phi_o, dchi_o); + cDpe = innerProduct(chi_e, dphi_e); + cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e compact " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o compact " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) compact " << pDce - conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) compact " << pDco - conj(cDpe) << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; std::cout << GridLogMessage << "= Test MeeInv Mee = 1 (if csw!=0) " << std::endl; std::cout << GridLogMessage << "==============================================================" << std::endl; @@ -169,7 +203,21 @@ int main(int argc, char **argv) setCheckerboard(phi, phi_o); err = phi - chi; - std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.Mooee(chi_e, src_e); + Dwc_compact.MooeeInv(src_e, phi_e); + + Dwc_compact.Mooee(chi_o, src_o); + Dwc_compact.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "==============================================================" << std::endl; std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1 (if csw!=0) " << std::endl; @@ -188,7 +236,21 @@ int main(int argc, char **argv) setCheckerboard(phi, phi_o); err = phi - chi; - std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.MooeeDag(chi_e, src_e); + Dwc_compact.MooeeInvDag(src_e, phi_e); + + Dwc_compact.MooeeDag(chi_o, src_o); + Dwc_compact.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "==============================================================" << std::endl; std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1 (if csw!=0) " << std::endl; @@ -207,7 +269,21 @@ int main(int argc, char **argv) setCheckerboard(phi, phi_o); err = phi - chi; - std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.MooeeDag(chi_e, src_e); + Dwc_compact.MooeeInv(src_e, phi_e); + + Dwc_compact.MooeeDag(chi_o, src_o); + Dwc_compact.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "================================================================" << std::endl; std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning " << std::endl; @@ -249,7 +325,7 @@ int main(int argc, char **argv) ///////////////// WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params); - Dwc_prime.ImportGauge(U_prime); + CompactWilsonCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params); tmp = Omega * src; pickCheckerboard(Even, src_e, tmp); @@ -262,7 +338,37 @@ int main(int argc, char **argv) setCheckerboard(phi, phi_o); err = chi - adj(Omega) * phi; - std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + tmp = Zero(); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc_compact.Mooee(src_e, chi_e); + Dwc_compact.Mooee(src_o, chi_o); + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + tmp = Omega * src; + pickCheckerboard(Even, src_e, tmp); + pickCheckerboard(Odd, src_o, tmp); + + Dwc_compact_prime.Mooee(src_e, phi_e); + Dwc_compact_prime.Mooee(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "=================================================================" << std::endl; std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning " << std::endl; @@ -272,7 +378,6 @@ int main(int argc, char **argv) phi = Zero(); WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params); - Dw.ImportGauge(Umu); Dw.M(src, result); Dwc.M(src, chi); @@ -280,13 +385,24 @@ int main(int argc, char **argv) Dwc_prime.M(Omega * src, phi); WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params); - Dw_prime.ImportGauge(U_prime); Dw_prime.M(Omega * src, result2); + err = result - adj(Omega) * result2; + std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); err = chi - adj(Omega) * phi; - err2 = result - adj(Omega) * result2; - std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl; - std::cout << GridLogMessage << "norm diff WilsonClover " << norm2(err2) << std::endl; + std::cout << GridLogMessage << "norm diff WilsonClover " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + + Dwc_compact.M(src, chi); + Dwc_compact_prime.M(Omega * src, phi); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff CompactWilsonClover " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "==========================================================" << std::endl; std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson " << std::endl; @@ -296,7 +412,6 @@ int main(int argc, char **argv) phi = Zero(); err = Zero(); WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0 - Dwc_csw0.ImportGauge(Umu); pickCheckerboard(Even, phi_e, phi); pickCheckerboard(Odd, phi_o, phi); @@ -316,7 +431,34 @@ int main(int argc, char **argv) setCheckerboard(src, src_o); err = chi - phi; - std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + err = Zero(); + CompactWilsonCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0 + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dw.Mooee(src_e, chi_e); + Dw.Mooee(src_o, chi_o); + Dwc_compact_csw0.Mooee(src_e, phi_e); + Dwc_compact_csw0.Mooee(src_o, phi_o); + + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + err = chi - phi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); std::cout << GridLogMessage << "==========================================================" << std::endl; std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec " << std::endl; @@ -348,9 +490,41 @@ int main(int argc, char **argv) setCheckerboard(phi, phi_o); err = ref - phi; - std::cout << GridLogMessage << "ref (unpreconditioned operator) diff :" << norm2(ref) << std::endl; - std::cout << GridLogMessage << "phi (EO decomposition) diff :" << norm2(phi) << std::endl; - std::cout << GridLogMessage << "norm diff :" << norm2(err) << std::endl; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff : " << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff : " << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff : " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + err = Zero(); + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o) + + Dwc_compact.M(src, ref); // Reference result from the unpreconditioned operator + + // EO matrix + Dwc_compact.Mooee(src_e, chi_e); + Dwc_compact.Mooee(src_o, chi_o); + Dwc_compact.Meooe(src_o, phi_e); + Dwc_compact.Meooe(src_e, phi_o); + + phi_o += chi_o; + phi_e += chi_e; + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = ref - phi; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff compact : " << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff compact : " << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff compact : " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); Grid_finalize(); } diff --git a/tests/core/Test_wilson_conserved_current.cc b/tests/core/Test_wilson_conserved_current.cc new file mode 100644 index 00000000..3ee1a271 --- /dev/null +++ b/tests/core/Test_wilson_conserved_current.cc @@ -0,0 +1,253 @@ +/************************************************************************************* +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_cayley_cg.cc + +Copyright (C) 2022 + +Author: Peter Boyle +Author: Fabian Joswig + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace std; +using namespace Grid; + + +template +void TestConserved(What & Dw, + LatticeGaugeField &Umu, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + GridParallelRNG *RNG4); + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT, + Gamma::Algebra::Gamma5 + }; + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout< seeds5({5,6,7,8}); + GridParallelRNG RNG4(UGrid); + std::vector seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); + if( argc > 1 && argv[1][0] != '-' ) + { + std::cout<::HotConfiguration(RNG4,Umu); + } + + typename WilsonCloverFermionR::ImplParams params; + WilsonAnisotropyCoefficients anis; + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + + std::cout<(Dw,Umu,UGrid,UrbGrid,&RNG4); + + std::cout<(Dwc,Umu,UGrid,UrbGrid,&RNG4); + + std::cout<(Dwcc,Umu,UGrid,UrbGrid,&RNG4); + + std::cout<(Dewc,Umu,UGrid,UrbGrid,&RNG4); + + std::cout<(Dewcc,Umu,UGrid,UrbGrid,&RNG4); + + Grid_finalize(); +} + + + +template +void TestConserved(Action & Dw, + LatticeGaugeField &Umu, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + GridParallelRNG *RNG4) +{ + LatticePropagator phys_src(UGrid); + LatticePropagator seqsrc(UGrid); + LatticePropagator prop4(UGrid); + LatticePropagator Vector_mu(UGrid); + LatticeComplex SV (UGrid); + LatticeComplex VV (UGrid); + LatticePropagator seqprop(UGrid); + + SpinColourMatrix kronecker; kronecker=1.0; + Coordinate coor({0,0,0,0}); + phys_src=Zero(); + pokeSite(kronecker,phys_src,coor); + + ConjugateGradient CG(1.0e-16,100000); + SchurRedBlackDiagTwoSolve schur(CG); + ZeroGuesser zpg; + for(int s=0;s(src4,phys_src,s,c); + + LatticeFermion result4(UGrid); result4=Zero(); + schur(Dw,src4,result4,zpg); + std::cout<(prop4,result4,s,c); + } + } + + auto curr = Current::Vector; + const int mu_J=0; + const int t_J=0; + + LatticeComplex ph (UGrid); ph=1.0; + + Dw.SeqConservedCurrent(prop4, + seqsrc, + phys_src, + curr, + mu_J, + t_J, + t_J,// whole lattice + ph); + + for(int s=0;s(src4,seqsrc,s,c); + + LatticeFermion result4(UGrid); result4=Zero(); + schur(Dw,src4,result4,zpg); + + FermToProp(seqprop,result4,s,c); + } + } + + Gamma g5(Gamma::Algebra::Gamma5); + Gamma gT(Gamma::Algebra::GammaT); + + std::vector sumSV; + std::vector sumVV; + + Dw.ContractConservedCurrent(prop4,prop4,Vector_mu,phys_src,Current::Vector,Tdir); + + SV = trace(Vector_mu); // Scalar-Vector conserved current + VV = trace(gT*Vector_mu); // (local) Vector-Vector conserved current + + // Spatial sum + sliceSum(SV,sumSV,Tdir); + sliceSum(VV,sumVV,Tdir); + + const int Nt{static_cast(sumSV.size())}; + + std::cout< check_buf; + + test_S = trace(qSite*g); + test_V = trace(qSite*g*Gamma::gmu[mu_J]); + + Dw.ContractConservedCurrent(prop4,prop4,cur,phys_src,curr,mu_J); + + c = trace(cur*g); + sliceSum(c, check_buf, Tp); + check_S = TensorRemove(check_buf[t_J]); + + auto gmu=Gamma::gmu[mu_J]; + c = trace(cur*g*gmu); + sliceSum(c, check_buf, Tp); + check_V = TensorRemove(check_buf[t_J]); + + + std::cout< + +using namespace std; +using namespace Grid; + +int main(int argc, char **argv) +{ + Grid_init(&argc, &argv); + + auto latt_size = GridDefaultLatt(); + auto simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd()); + auto mpi_layout = GridDefaultMpi(); + GridCartesian Grid(latt_size, simd_layout, mpi_layout); + GridRedBlackCartesian RBGrid(&Grid); + + int threads = GridThread::GetThreads(); + std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REALF" << sizeof(RealF) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REALD" << sizeof(RealD) << std::endl; + std::cout << GridLogMessage << "Grid floating point word size is REAL" << sizeof(Real) << std::endl; + + std::vector seeds({1, 2, 3, 4}); + GridParallelRNG pRNG(&Grid); + pRNG.SeedFixedIntegers(seeds); + // pRNG.SeedFixedIntegers(std::vector({45,12,81,9}); + + typedef typename WilsonExpCloverFermionR::FermionField FermionField; + typename WilsonExpCloverFermionR::ImplParams params; + WilsonAnisotropyCoefficients anis; + + FermionField src(&Grid); + random(pRNG, src); + FermionField result(&Grid); + result = Zero(); + FermionField result2(&Grid); + result2 = Zero(); + FermionField ref(&Grid); + ref = Zero(); + FermionField tmp(&Grid); + tmp = Zero(); + FermionField err(&Grid); + err = Zero(); + FermionField phi(&Grid); + random(pRNG, phi); + FermionField chi(&Grid); + random(pRNG, chi); + LatticeGaugeField Umu(&Grid); + SU::HotConfiguration(pRNG, Umu); + std::vector U(4, &Grid); + + double tolerance = 1e-4; + + double volume = 1; + for (int mu = 0; mu < Nd; mu++) + { + volume = volume * latt_size[mu]; + } + + RealD mass = 0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + + WilsonExpCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + CompactWilsonExpCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params); + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + FermionField src_e(&RBGrid); + FermionField src_o(&RBGrid); + FermionField r_e(&RBGrid); + FermionField r_o(&RBGrid); + FermionField r_eo(&Grid); + pickCheckerboard(Even, src_e, src); + pickCheckerboard(Odd, src_o, src); + + Dwc.Meooe(src_e, r_o); + std::cout << GridLogMessage << "Applied Meo" << std::endl; + Dwc.Meooe(src_o, r_e); + std::cout << GridLogMessage << "Applied Moe" << std::endl; + Dwc.Dhop(src, ref, DaggerNo); + + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl; + assert(fabs(norm2(err)) < tolerance); + + + + Dwc_compact.Meooe(src_e, r_o); + std::cout << GridLogMessage << "Applied Meo" << std::endl; + Dwc_compact.Meooe(src_o, r_e); + std::cout << GridLogMessage << "Applied Moe" << std::endl; + Dwc_compact.Dhop(src, ref, DaggerNo); + + setCheckerboard(r_eo, r_o); + setCheckerboard(r_eo, r_e); + + err = ref - r_eo; + std::cout << GridLogMessage << "EO norm diff compact\t" << norm2(err) << " (" << norm2(ref) << " - " << norm2(r_eo) << ")" << std::endl; + assert(fabs(norm2(err)) < tolerance); + + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test Ddagger is the dagger of D by requiring " << std::endl; + std::cout << GridLogMessage << "= < phi | Deo | chi > * = < chi | Deo^dag| phi> " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + FermionField chi_e(&RBGrid); + FermionField chi_o(&RBGrid); + + FermionField dchi_e(&RBGrid); + FermionField dchi_o(&RBGrid); + + FermionField phi_e(&RBGrid); + FermionField phi_o(&RBGrid); + + FermionField dphi_e(&RBGrid); + FermionField dphi_o(&RBGrid); + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Meooe(chi_e, dchi_o); + Dwc.Meooe(chi_o, dchi_e); + Dwc.MeooeDag(phi_e, dphi_o); + Dwc.MeooeDag(phi_o, dphi_e); + + ComplexD pDce = innerProduct(phi_e, dchi_e); + ComplexD pDco = innerProduct(phi_o, dchi_o); + ComplexD cDpe = innerProduct(chi_e, dphi_e); + ComplexD cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) " << pDce - conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) " << pDco - conj(cDpe) << std::endl; + + Dwc_compact.Meooe(chi_e, dchi_o); + Dwc_compact.Meooe(chi_o, dchi_e); + Dwc_compact.MeooeDag(phi_e, dphi_o); + Dwc_compact.MeooeDag(phi_o, dphi_e); + + pDce = innerProduct(phi_e, dchi_e); + pDco = innerProduct(phi_o, dchi_o); + cDpe = innerProduct(chi_e, dphi_e); + cDpo = innerProduct(chi_o, dphi_o); + + std::cout << GridLogMessage << "e compact " << pDce << " " << cDpe << std::endl; + std::cout << GridLogMessage << "o compact " << pDco << " " << cDpo << std::endl; + + std::cout << GridLogMessage << "pDce - conj(cDpo) compact " << pDce - conj(cDpo) << std::endl; + std::cout << GridLogMessage << "pDco - conj(cDpe) compact " << pDco - conj(cDpe) << std::endl; + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv Mee = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.Mooee(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.Mooee(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.Mooee(chi_e, src_e); + Dwc_compact.MooeeInv(src_e, phi_e); + + Dwc_compact.Mooee(chi_o, src_o); + Dwc_compact.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeDag MeeInvDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInvDag(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.MooeeDag(chi_e, src_e); + Dwc_compact.MooeeInvDag(src_e, phi_e); + + Dwc_compact.MooeeDag(chi_o, src_o); + Dwc_compact.MooeeInvDag(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "==============================================================" << std::endl; + std::cout << GridLogMessage << "= Test MeeInv MeeDag = 1 (if csw!=0) " << std::endl; + std::cout << GridLogMessage << "==============================================================" << std::endl; + + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dwc.MooeeDag(chi_e, src_e); + Dwc.MooeeInv(src_e, phi_e); + + Dwc.MooeeDag(chi_o, src_o); + Dwc.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Dwc_compact.MooeeDag(chi_e, src_e); + Dwc_compact.MooeeInv(src_e, phi_e); + + Dwc_compact.MooeeDag(chi_o, src_o); + Dwc_compact.MooeeInv(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = phi - chi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term with EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = Zero(); + phi = Zero(); + tmp = Zero(); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + ////////////////////// Gauge Transformation + std::vector seeds2({5, 6, 7, 8}); + GridParallelRNG pRNG2(&Grid); + pRNG2.SeedFixedIntegers(seeds2); + LatticeColourMatrix Omega(&Grid); + LatticeColourMatrix ShiftedOmega(&Grid); + LatticeGaugeField U_prime(&Grid); + U_prime = Zero(); + LatticeColourMatrix U_prime_mu(&Grid); + U_prime_mu = Zero(); + SU::LieRandomize(pRNG2, Omega, 1.0); + for (int mu = 0; mu < Nd; mu++) + { + U[mu] = peekLorentz(Umu, mu); + ShiftedOmega = Cshift(Omega, mu, 1); + U_prime_mu = Omega * U[mu] * adj(ShiftedOmega); + pokeLorentz(U_prime, U_prime_mu, mu); + } + ///////////////// + + WilsonExpCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params); + CompactWilsonExpCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params); + + tmp = Omega * src; + pickCheckerboard(Even, src_e, tmp); + pickCheckerboard(Odd, src_o, tmp); + + Dwc_prime.Mooee(src_e, phi_e); + Dwc_prime.Mooee(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + tmp = Zero(); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + + Dwc_compact.Mooee(src_e, chi_e); + Dwc_compact.Mooee(src_o, chi_o); + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + tmp = Omega * src; + pickCheckerboard(Even, src_e, tmp); + pickCheckerboard(Odd, src_o, tmp); + + Dwc_compact_prime.Mooee(src_e, phi_e); + Dwc_compact_prime.Mooee(src_o, phi_o); + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "=================================================================" << std::endl; + std::cout << GridLogMessage << "= Testing gauge covariance Clover term w/o EO preconditioning " << std::endl; + std::cout << GridLogMessage << "================================================================" << std::endl; + + chi = Zero(); + phi = Zero(); + + WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params); + + Dw.M(src, result); + Dwc.M(src, chi); + + Dwc_prime.M(Omega * src, phi); + + WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params); + Dw_prime.M(Omega * src, result2); + + err = result - adj(Omega) * result2; + std::cout << GridLogMessage << "norm diff Wilson " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff WilsonExpClover " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + + Dwc_compact.M(src, chi); + Dwc_compact_prime.M(Omega * src, phi); + + err = chi - adj(Omega) * phi; + std::cout << GridLogMessage << "norm diff CompactWilsonExpClover " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing Mooee(csw=0) Clover to reproduce Mooee Wilson " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = Zero(); + phi = Zero(); + err = Zero(); + WilsonExpCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0 + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dw.Mooee(src_e, chi_e); + Dw.Mooee(src_o, chi_o); + Dwc_csw0.Mooee(src_e, phi_e); + Dwc_csw0.Mooee(src_o, phi_o); + + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + err = chi - phi; + std::cout << GridLogMessage << "norm diff " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + err = Zero(); + CompactWilsonExpCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0 + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + Dw.Mooee(src_e, chi_e); + Dw.Mooee(src_o, chi_o); + Dwc_compact_csw0.Mooee(src_e, phi_e); + Dwc_compact_csw0.Mooee(src_o, phi_o); + + setCheckerboard(chi, chi_e); + setCheckerboard(chi, chi_o); + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + setCheckerboard(src, src_e); + setCheckerboard(src, src_o); + + err = chi - phi; + std::cout << GridLogMessage << "norm diff compact " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + std::cout << GridLogMessage << "==========================================================" << std::endl; + std::cout << GridLogMessage << "= Testing EO operator is equal to the unprec " << std::endl; + std::cout << GridLogMessage << "==========================================================" << std::endl; + + chi = Zero(); + phi = Zero(); + err = Zero(); + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o) + + Dwc.M(src, ref); // Reference result from the unpreconditioned operator + + // EO matrix + Dwc.Mooee(src_e, chi_e); + Dwc.Mooee(src_o, chi_o); + Dwc.Meooe(src_o, phi_e); + Dwc.Meooe(src_e, phi_o); + + phi_o += chi_o; + phi_e += chi_e; + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = ref - phi; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff : " << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff : " << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff : " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + chi = Zero(); + phi = Zero(); + err = Zero(); + + pickCheckerboard(Even, phi_e, phi); + pickCheckerboard(Odd, phi_o, phi); + pickCheckerboard(Even, chi_e, chi); + pickCheckerboard(Odd, chi_o, chi); + + // M phi = (Mooee src_e + Meooe src_o , Meooe src_e + Mooee src_o) + + Dwc_compact.M(src, ref); // Reference result from the unpreconditioned operator + + // EO matrix + Dwc_compact.Mooee(src_e, chi_e); + Dwc_compact.Mooee(src_o, chi_o); + Dwc_compact.Meooe(src_o, phi_e); + Dwc_compact.Meooe(src_e, phi_o); + + phi_o += chi_o; + phi_e += chi_e; + + setCheckerboard(phi, phi_e); + setCheckerboard(phi, phi_o); + + err = ref - phi; + std::cout << GridLogMessage << "ref (unpreconditioned operator) diff compact : " << norm2(ref) << std::endl; + std::cout << GridLogMessage << "phi (EO decomposition) diff compact : " << norm2(phi) << std::endl; + std::cout << GridLogMessage << "norm diff compact : " << norm2(err) << std::endl; + assert(fabs(norm2(err)) < tolerance); + + Grid_finalize(); +} diff --git a/tests/hmc/Test_hmc_ScalarActionNxN.cc b/tests/hmc/Test_hmc_ScalarActionNxN.cc index 0c398306..726ecd4a 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.LoadScidacCheckpointer(CPparams, SPar); this breaks for compilation without lime 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..8cd74791 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); + // Store metadata in the Scidac checkpointer - obviously breaks without LIME + //TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar); RNGModuleParameters RNGpar(Reader); TheHMC.Resources.SetRNGSeeds(RNGpar); diff --git a/tests/lanczos/Test_compressed_lanczos.cc b/tests/lanczos/Test_compressed_lanczos.cc index 3fa3c490..d7d0d52d 100644 --- a/tests/lanczos/Test_compressed_lanczos.cc +++ b/tests/lanczos/Test_compressed_lanczos.cc @@ -37,6 +37,8 @@ Author: Peter Boyle using namespace std; using namespace Grid; +#ifdef HAVE_LIME + template class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos { @@ -249,3 +251,11 @@ int main (int argc, char ** argv) { Grid_finalize(); } +#else + +int main( void ) +{ + return 0 ; +} + +#endif // HAVE_LIME_H diff --git a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc index ecd0e629..669a7b6d 100644 --- a/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc +++ b/tests/lanczos/Test_dwf_compressed_lanczos_reorg.cc @@ -36,7 +36,8 @@ Author: Peter Boyle using namespace std; using namespace Grid; - ; + +#ifdef HAVE_LIME template class LocalCoherenceLanczosScidac : public LocalCoherenceLanczos @@ -250,3 +251,11 @@ int main (int argc, char ** argv) { Grid_finalize(); } +#else + +int main( void ) +{ + return 0 ; +} + +#endif diff --git a/tests/solver/Test_coarse_even_odd.cc b/tests/solver/Test_coarse_even_odd.cc index dfbab747..c7127121 100644 --- a/tests/solver/Test_coarse_even_odd.cc +++ b/tests/solver/Test_coarse_even_odd.cc @@ -93,8 +93,16 @@ int main(int argc, char** argv) { // Setup of Dirac Matrix and Operator // ///////////////////////////////////////////////////////////////////////////// - LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu); - + LatticeGaugeField Umu(Grid_f); +#if (Nc==2) + SU2::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==3) + SU3::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==4) + SU4::HotConfiguration(pRNG_f, Umu); +#elif (defined Nc==5) + SU5::HotConfiguration(pRNG_f, Umu); +#endif RealD checkTolerance = (getPrecision::value == 1) ? 1e-7 : 1e-15; RealD mass = -0.30; diff --git a/tests/solver/Test_dwf_mrhs_cg.cc b/tests/solver/Test_dwf_mrhs_cg.cc index b912ba4f..4242f64b 100644 --- a/tests/solver/Test_dwf_mrhs_cg.cc +++ b/tests/solver/Test_dwf_mrhs_cg.cc @@ -34,6 +34,7 @@ using namespace Grid; int main (int argc, char ** argv) { +#ifdef HAVE_LIME typedef typename DomainWallFermionR::FermionField FermionField; typedef typename DomainWallFermionR::ComplexField ComplexField; typename DomainWallFermionR::ImplParams params; @@ -237,4 +238,5 @@ int main (int argc, char ** argv) } Grid_finalize(); +#endif // HAVE_LIME } diff --git a/tests/solver/Test_wilsonclover_cg_prec.cc b/tests/solver/Test_wilsonclover_cg_prec.cc index 5ef2dc09..abf64a1f 100644 --- a/tests/solver/Test_wilsonclover_cg_prec.cc +++ b/tests/solver/Test_wilsonclover_cg_prec.cc @@ -71,7 +71,12 @@ int main (int argc, char ** argv) RealD mass = -0.1; RealD csw_r = 1.0; RealD csw_t = 1.0; + RealD cF = 1.0; WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); + WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t); + CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); + // HermitianOperator HermOp(Dw); // ConjugateGradient CG(1.0e-8,10000); @@ -80,12 +85,28 @@ int main (int argc, char ** argv) LatticeFermion src_o(&RBGrid); LatticeFermion result_o(&RBGrid); pickCheckerboard(Odd,src_o,src); - result_o=Zero(); - SchurDiagMooeeOperator HermOpEO(Dw); ConjugateGradient CG(1.0e-8,10000); + + std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl; + SchurDiagMooeeOperator HermOpEO(Dw); + result_o=Zero(); CG(HermOpEO,src_o,result_o); - + + std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl; + SchurDiagMooeeOperator HermOpEO_compact(Dw_compact); + result_o=Zero(); + CG(HermOpEO_compact,src_o,result_o); + + std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl; + SchurDiagMooeeOperator HermOpEO_exp(Dwe); + result_o=Zero(); + CG(HermOpEO_exp,src_o,result_o); + + std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl; + SchurDiagMooeeOperator HermOpEO_exp_compact(Dwe_compact); + result_o=Zero(); + CG(HermOpEO_exp_compact,src_o,result_o); Grid_finalize(); } diff --git a/tests/solver/Test_wilsonclover_cg_schur.cc b/tests/solver/Test_wilsonclover_cg_schur.cc index 567a8283..50d06af7 100644 --- a/tests/solver/Test_wilsonclover_cg_schur.cc +++ b/tests/solver/Test_wilsonclover_cg_schur.cc @@ -60,18 +60,36 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(&Grid); SU::HotConfiguration(pRNG,Umu); LatticeFermion src(&Grid); random(pRNG,src); - LatticeFermion result(&Grid); result=Zero(); - LatticeFermion resid(&Grid); - - RealD mass = -0.1; - RealD csw_r = 1.0; - RealD csw_t = 1.0; - WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + LatticeFermion result(&Grid); + LatticeFermion resid(&Grid); ConjugateGradient CG(1.0e-8,10000); SchurRedBlackDiagMooeeSolve SchurSolver(CG); + RealD mass = -0.1; + RealD csw_r = 1.0; + RealD csw_t = 1.0; + RealD cF = 1.0; + + std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl; + WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + result=Zero(); SchurSolver(Dw,src,result); + + std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl; + CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); + result=Zero(); + SchurSolver(Dw_compact,src,result); + + std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl; + WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t); + result=Zero(); + SchurSolver(Dwe,src,result); + + std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl; + CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); + result=Zero(); + SchurSolver(Dwe_compact,src,result); Grid_finalize(); } diff --git a/tests/solver/Test_wilsonclover_cg_unprec.cc b/tests/solver/Test_wilsonclover_cg_unprec.cc index 755d80e1..2a859f11 100644 --- a/tests/solver/Test_wilsonclover_cg_unprec.cc +++ b/tests/solver/Test_wilsonclover_cg_unprec.cc @@ -59,7 +59,7 @@ int main (int argc, char ** argv) LatticeFermion src(&Grid); random(pRNG,src); RealD nrm = norm2(src); - LatticeFermion result(&Grid); result=Zero(); + LatticeFermion result(&Grid); LatticeGaugeField Umu(&Grid); SU::HotConfiguration(pRNG,Umu); double volume=1; @@ -70,11 +70,34 @@ int main (int argc, char ** argv) RealD mass = -0.1; RealD csw_r = 1.0; RealD csw_t = 1.0; + RealD cF = 1.0; WilsonCloverFermionR Dw(Umu, Grid, RBGrid, mass, csw_r, csw_t); + CompactWilsonCloverFermionR Dw_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); + WilsonExpCloverFermionR Dwe(Umu, Grid, RBGrid, mass, csw_r, csw_t); + CompactWilsonExpCloverFermionR Dwe_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 0.0); - MdagMLinearOperator HermOp(Dw); ConjugateGradient CG(1.0e-8,10000); + + std::cout << GridLogMessage << "Testing Wilson Clover" << std::endl; + MdagMLinearOperator HermOp(Dw); + result=Zero(); CG(HermOp,src,result); + std::cout << GridLogMessage << "Testing Compact Wilson Clover" << std::endl; + MdagMLinearOperator HermOp_compact(Dw_compact); + result=Zero(); + CG(HermOp_compact,src,result); + + + std::cout << GridLogMessage << "Testing Wilson Exp Clover" << std::endl; + MdagMLinearOperator HermOp_exp(Dwe); + result=Zero(); + CG(HermOp_exp,src,result); + + std::cout << GridLogMessage << "Testing Compact Wilson Exp Clover" << std::endl; + MdagMLinearOperator HermOp_exp_compact(Dwe_compact); + result=Zero(); + CG(HermOp_exp_compact,src,result); + Grid_finalize(); }