portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 03:05:55 +01:00
Code Releases Activity

Charge conjugation boundary conditions for gauge fields implemented as a policy

Browse Source 
class, changing the nature of covariant Cshifts used in
plaquettes, rectangles and staples.

As a result same code is used for the plaq and rect action independent of the BC type.

Should probably isolate the BC in a separate class that Gimpl takes as a template param.
Do the same with smearing policies.

This would then allow composition of BC with smearing etc....
This commit is contained in:
paboyle 2016-01-02 13:37:25 +00:00
parent 145a295231
commit 5a80930dd2
20 changed files with 672 additions and 213 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
lib/Make.inc
View File

@ -1,4 +1,4 @@
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/NerscCheckpointer.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/GaugeImpl.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/NerscCheckpointer.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./PerfCount.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsAsm.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonTMFermion.cc ./qcd/hmc/HMC.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./serialisation/BinaryIO.cc ./serialisation/TextIO.cc ./serialisation/XmlIO.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./PerfCount.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsAsm.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonTMFermion.cc ./qcd/hmc/HMC.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./serialisation/BinaryIO.cc ./serialisation/TextIO.cc ./serialisation/XmlIO.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

6
lib/parallelIO/NerscIO.h
View File

@ -51,7 +51,6 @@ class NerscField {
std::string floating_point; std::string floating_point;
}; };
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
// Bit and Physical Checksumming and QA of data // Bit and Physical Checksumming and QA of data
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
@ -69,8 +68,9 @@ inline void NerscGrid(GridBase *grid,NerscField &header)
template<class GaugeField> template<class GaugeField>
inline void NerscStatistics(GaugeField & data,NerscField &header) inline void NerscStatistics(GaugeField & data,NerscField &header)
{ {
header.link_trace=Grid::QCD::WilsonLoops<GaugeField>::linkTrace(data); // How to convert data precision etc...
header.plaquette =Grid::QCD::WilsonLoops<GaugeField>::avgPlaquette(data); header.link_trace=Grid::QCD::WilsonLoops<PeriodicGimplR>::linkTrace(data);
header.plaquette =Grid::QCD::WilsonLoops<PeriodicGimplR>::avgPlaquette(data);
} }
inline void NerscMachineCharacteristics(NerscField &header) inline void NerscMachineCharacteristics(NerscField &header)

1
lib/qcd/QCD.h
View File

@ -425,7 +425,6 @@ namespace QCD {
#include <qcd/spin/TwoSpinor.h> #include <qcd/spin/TwoSpinor.h>
#include <qcd/utils/LinalgUtils.h> #include <qcd/utils/LinalgUtils.h>
#include <qcd/utils/CovariantCshift.h> #include <qcd/utils/CovariantCshift.h>
#include <qcd/utils/WilsonLoops.h>
#include <qcd/utils/SUn.h> #include <qcd/utils/SUn.h>
#include <qcd/action/Actions.h> #include <qcd/action/Actions.h>
#include <qcd/hmc/integrators/Integrator.h> #include <qcd/hmc/integrators/Integrator.h>

54
lib/qcd/action/Actions.h
View File

@ -1,11 +1,6 @@
#ifndef GRID_QCD_ACTIONS_H #ifndef GRID_QCD_ACTIONS_H
#define GRID_QCD_ACTIONS_H #define GRID_QCD_ACTIONS_H
// Some reorganisation likely required as both Chroma and IroIro
// are separating the concept of the operator from that of action.
//
// The FermAction contains methods to create
// * Linear operators (Hermitian and non-hermitian) .. my LinearOperator // * Linear operators (Hermitian and non-hermitian) .. my LinearOperator
// * System solvers (Hermitian and non-hermitian) .. my OperatorFunction // * System solvers (Hermitian and non-hermitian) .. my OperatorFunction
// * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction // * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction
@ -19,6 +14,9 @@
//////////////////////////////////////////// ////////////////////////////////////////////
// Utility functions // Utility functions
//////////////////////////////////////////// ////////////////////////////////////////////
#include <qcd/action/gauge/GaugeImpl.h>
#include <qcd/utils/WilsonLoops.h>
#include <qcd/action/fermion/WilsonCompressor.h> //used by all wilson type fermions #include <qcd/action/fermion/WilsonCompressor.h> //used by all wilson type fermions
#include <qcd/action/fermion/FermionOperatorImpl.h> #include <qcd/action/fermion/FermionOperatorImpl.h>
#include <qcd/action/fermion/FermionOperator.h> #include <qcd/action/fermion/FermionOperator.h>
@ -29,20 +27,37 @@
//////////////////////////////////////////// ////////////////////////////////////////////
#include <qcd/action/gauge/WilsonGaugeAction.h> #include <qcd/action/gauge/WilsonGaugeAction.h>
#include <qcd/action/gauge/PlaqPlusRectangleAction.h> #include <qcd/action/gauge/PlaqPlusRectangleAction.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
typedef WilsonGaugeAction<LatticeGaugeField> WilsonGaugeActionR;
typedef WilsonGaugeAction<LatticeGaugeFieldF> WilsonGaugeActionF; typedef WilsonGaugeAction<PeriodicGimplR> WilsonGaugeActionR;
typedef WilsonGaugeAction<LatticeGaugeFieldD> WilsonGaugeActionD; typedef WilsonGaugeAction<PeriodicGimplF> WilsonGaugeActionF;
typedef PlaqPlusRectangleAction<LatticeGaugeField> PlaqPlusRectangleActionR; typedef WilsonGaugeAction<PeriodicGimplD> WilsonGaugeActionD;
typedef PlaqPlusRectangleAction<LatticeGaugeFieldF> PlaqPlusRectangleActionF; typedef PlaqPlusRectangleAction<PeriodicGimplR> PlaqPlusRectangleActionR;
typedef PlaqPlusRectangleAction<LatticeGaugeFieldD> PlaqPlusRectangleActionD; typedef PlaqPlusRectangleAction<PeriodicGimplF> PlaqPlusRectangleActionF;
typedef IwasakiGaugeAction<LatticeGaugeField> IwasakiGaugeActionR; typedef PlaqPlusRectangleAction<PeriodicGimplD> PlaqPlusRectangleActionD;
typedef IwasakiGaugeAction<LatticeGaugeFieldF> IwasakiGaugeActionF; typedef IwasakiGaugeAction<PeriodicGimplR> IwasakiGaugeActionR;
typedef IwasakiGaugeAction<LatticeGaugeFieldD> IwasakiGaugeActionD; typedef IwasakiGaugeAction<PeriodicGimplF> IwasakiGaugeActionF;
typedef SymanzikGaugeAction<LatticeGaugeField> SymanzikGaugeActionR; typedef IwasakiGaugeAction<PeriodicGimplD> IwasakiGaugeActionD;
typedef SymanzikGaugeAction<LatticeGaugeFieldF> SymanzikGaugeActionF; typedef SymanzikGaugeAction<PeriodicGimplR> SymanzikGaugeActionR;
typedef SymanzikGaugeAction<LatticeGaugeFieldD> SymanzikGaugeActionD; typedef SymanzikGaugeAction<PeriodicGimplF> SymanzikGaugeActionF;
typedef SymanzikGaugeAction<PeriodicGimplD> SymanzikGaugeActionD;
typedef WilsonGaugeAction<ConjugateGimplR> ConjugateWilsonGaugeActionR;
typedef WilsonGaugeAction<ConjugateGimplF> ConjugateWilsonGaugeActionF;
typedef WilsonGaugeAction<ConjugateGimplD> ConjugateWilsonGaugeActionD;
typedef PlaqPlusRectangleAction<ConjugateGimplR> ConjugatePlaqPlusRectangleActionR;
typedef PlaqPlusRectangleAction<ConjugateGimplF> ConjugatePlaqPlusRectangleActionF;
typedef PlaqPlusRectangleAction<ConjugateGimplD> ConjugatePlaqPlusRectangleActionD;
typedef IwasakiGaugeAction<ConjugateGimplR> ConjugateIwasakiGaugeActionR;
typedef IwasakiGaugeAction<ConjugateGimplF> ConjugateIwasakiGaugeActionF;
typedef IwasakiGaugeAction<ConjugateGimplD> ConjugateIwasakiGaugeActionD;
typedef SymanzikGaugeAction<ConjugateGimplR> ConjugateSymanzikGaugeActionR;
typedef SymanzikGaugeAction<ConjugateGimplF> ConjugateSymanzikGaugeActionF;
typedef SymanzikGaugeAction<ConjugateGimplD> ConjugateSymanzikGaugeActionD;
}} }}
//////////////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////////
@ -177,11 +192,6 @@ typedef DomainWallFermion<GparityWilsonImplD> GparityDomainWallFermionD;
#include <qcd/action/pseudofermion/TwoFlavourEvenOdd.h> #include <qcd/action/pseudofermion/TwoFlavourEvenOdd.h>
#include <qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h> #include <qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h>
//IroIro inserted general "Nf" param; could also be done,
//but not clear why unless into large Nf BSM studies
//Even there, don't want the explicit (2) on power denominator
//if even number of flavours, so further generalised interface
//would be required but easy.
#include <qcd/action/pseudofermion/OneFlavourRational.h> #include <qcd/action/pseudofermion/OneFlavourRational.h>
#include <qcd/action/pseudofermion/OneFlavourRationalRatio.h> #include <qcd/action/pseudofermion/OneFlavourRationalRatio.h>
#include <qcd/action/pseudofermion/OneFlavourEvenOddRational.h> #include <qcd/action/pseudofermion/OneFlavourEvenOddRational.h>

44
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -58,38 +58,6 @@ namespace Grid {
// } // }
////////////////////////////////////////////// //////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
// Implementation dependent gauge types
////////////////////////////////////////////////////////////////////////
#define INHERIT_IMPL_TYPES(Base) \
INHERIT_GIMPL_TYPES(Base)\
INHERIT_FIMPL_TYPES(Base)
#define INHERIT_GIMPL_TYPES(GImpl) \
typedef typename GImpl::Simd Simd;\
typedef typename GImpl::GaugeLinkField GaugeLinkField;\
typedef typename GImpl::GaugeField GaugeField;
// Composition with smeared link, bc's etc.. probably need multiple inheritance
// Variable precision "S" and variable Nc
template<class S,int Nrepresentation=Nc>
class ImplGauge {
public:
typedef S Simd;
template<typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
template<typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd >;
typedef iImplGaugeLink <Simd> SiteGaugeLink;
typedef iImplGaugeField <Simd> SiteGaugeField;
typedef Lattice<SiteGaugeLink> GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
typedef Lattice<SiteGaugeField> GaugeField;
};
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Implementation dependent fermion types // Implementation dependent fermion types
@ -104,14 +72,18 @@ namespace Grid {
typedef typename Impl::StencilImpl StencilImpl; \ typedef typename Impl::StencilImpl StencilImpl; \
typedef typename Impl::ImplParams ImplParams; typedef typename Impl::ImplParams ImplParams;
#define INHERIT_IMPL_TYPES(Base) \
INHERIT_GIMPL_TYPES(Base)\
INHERIT_FIMPL_TYPES(Base)
/////// ///////
// Single flavour four spinors with colour index // Single flavour four spinors with colour index
/////// ///////
template<class S,int Nrepresentation=Nc> template<class S,int Nrepresentation=Nc>
class WilsonImpl : public ImplGauge<S,Nrepresentation> { class WilsonImpl : public PeriodicGaugeImpl<S,Nrepresentation> {
public: public:
typedef ImplGauge<S,Nrepresentation> Gimpl; typedef PeriodicGaugeImpl<S,Nrepresentation> Gimpl;
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);
@ -180,10 +152,10 @@ PARALLEL_FOR_LOOP
//////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////
template<class S,int Nrepresentation> template<class S,int Nrepresentation>
class GparityWilsonImpl : public ImplGauge<S,Nrepresentation> { class GparityWilsonImpl : public ConjugateGaugeImpl<S,Nrepresentation> {
public: public:
typedef ImplGauge<S,Nrepresentation> Gimpl; typedef ConjugateGaugeImpl<S,Nrepresentation> Gimpl;
INHERIT_GIMPL_TYPES(Gimpl); INHERIT_GIMPL_TYPES(Gimpl);

150
lib/qcd/action/gauge/GaugeImpl.h Normal file
View File

@ -0,0 +1,150 @@
#ifndef GRID_QCD_GAUGE_IMPL_H
#define GRID_QCD_GAUGE_IMPL_H
namespace Grid {
namespace QCD {
////////////////////////////////////////////////////////////////////////
// Implementation dependent gauge types
////////////////////////////////////////////////////////////////////////
template<class Gimpl> class WilsonLoops;
#define INHERIT_GIMPL_TYPES(GImpl) \
typedef typename GImpl::Simd Simd;\
typedef typename GImpl::GaugeLinkField GaugeLinkField;\
typedef typename GImpl::GaugeField GaugeField;
// Composition with smeared link, bc's etc.. probably need multiple inheritance
// Variable precision "S" and variable Nc
template<class S,int Nrepresentation=Nc>
class PeriodicGaugeImpl {
public:
typedef S Simd;
template<typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
template<typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd >;
typedef iImplGaugeLink <Simd> SiteGaugeLink;
typedef iImplGaugeField <Simd> SiteGaugeField;
typedef Lattice<SiteGaugeLink> GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
typedef Lattice<SiteGaugeField> GaugeField;
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Support needed for the assembly of loops including all boundary condition effects such as conjugate bcs
////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class covariant> static inline
Lattice<covariant> CovShiftForward (const GaugeLinkField &Link, int mu, const Lattice<covariant> &field) {
return PeriodicBC::CovShiftForward(Link,mu,field);
}
template<class covariant> static inline
Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,const Lattice<covariant> &field) {
return PeriodicBC::CovShiftBackward(Link,mu,field);
}
static inline
GaugeLinkField CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
return Cshift(adj(Link),mu,-1);
}
static inline
GaugeLinkField CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {
return Link;
}
static inline
GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
return Cshift(Link,mu,1);
}
static inline bool isPeriodicGaugeField(void) {
return true;
}
};
// Composition with smeared link, bc's etc.. probably need multiple inheritance
// Variable precision "S" and variable Nc
template<class S,int Nrepresentation=Nc>
class ConjugateGaugeImpl {
public:
typedef S Simd;
template<typename vtype> using iImplGaugeLink = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
template<typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd >;
typedef iImplGaugeLink <Simd> SiteGaugeLink;
typedef iImplGaugeField <Simd> SiteGaugeField;
typedef Lattice<SiteGaugeLink> GaugeLinkField; // bit ugly naming; polarised gauge field, lorentz... all ugly
typedef Lattice<SiteGaugeField> GaugeField;
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Support needed for the assembly of loops including all boundary condition effects such as Gparity.
////////////////////////////////////////////////////////////////////////////////////////////////////////////
template<class covariant> static
Lattice<covariant> CovShiftForward (const GaugeLinkField &Link, int mu, const Lattice<covariant> &field) {
return GparityBC::CovShiftForward(Link,mu,field);
}
template<class covariant> static
Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,const Lattice<covariant> &field) {
return GparityBC::CovShiftBackward(Link,mu,field);
}
static inline
GaugeLinkField CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) {
GridBase *grid = Link._grid;
int Lmu = grid->GlobalDimensions()[mu]-1;
Lattice<iScalar<vInteger> > coor(grid); LatticeCoordinate(coor,mu);
GaugeLinkField tmp (grid);
tmp=adj(Link);
tmp = where(coor==Lmu,conjugate(tmp),tmp);
return Cshift(tmp,mu,-1);// moves towards positive mu
}
static inline
GaugeLinkField CovShiftIdentityForward(const GaugeLinkField &Link, int mu) {
return Link;
}
static inline
GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) {
GridBase *grid = Link._grid;
int Lmu = grid->GlobalDimensions()[mu]-1;
Lattice<iScalar<vInteger> > coor(grid); LatticeCoordinate(coor,mu);
GaugeLinkField tmp (grid);
tmp=Cshift(Link,mu,1);
tmp=where(coor==Lmu,conjugate(tmp),tmp);
return tmp;
}
static inline bool isPeriodicGaugeField(void) {
return false;
}
};
typedef PeriodicGaugeImpl<vComplex ,Nc> PeriodicGimplR; // Real.. whichever prec
typedef PeriodicGaugeImpl<vComplexF,Nc> PeriodicGimplF; // Float
typedef PeriodicGaugeImpl<vComplexD,Nc> PeriodicGimplD; // Double
typedef ConjugateGaugeImpl<vComplex ,Nc> ConjugateGimplR; // Real.. whichever prec
typedef ConjugateGaugeImpl<vComplexF,Nc> ConjugateGimplF; // Float
typedef ConjugateGaugeImpl<vComplexD,Nc> ConjugateGimplD; // Double
}
}
#endif

46
lib/qcd/action/gauge/PlaqPlusRectangleAction.h
View File

@ -7,11 +7,11 @@ namespace Grid{
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// PlaqPlusRectangleActoin // PlaqPlusRectangleActoin
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
template<class GaugeField> template<class Gimpl>
class PlaqPlusRectangleAction : public Action<GaugeField> { class PlaqPlusRectangleAction : public Action<typename Gimpl::GaugeField> {
public: public:
typedef LorentzScalar<GaugeField> GaugeLinkField; INHERIT_GIMPL_TYPES(Gimpl);
private: private:
RealD c_plaq; RealD c_plaq;
@ -25,8 +25,8 @@ namespace Grid{
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
RealD vol = U._grid->gSites(); RealD vol = U._grid->gSites();
RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U); RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
RealD rect = WilsonLoops<GaugeField>::avgRectangle(U); RealD rect = WilsonLoops<Gimpl>::avgRectangle(U);
RealD action=c_plaq*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5 RealD action=c_plaq*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5
+c_rect*(1.0 -rect)*(Nd*(Nd-1.0))*vol; +c_rect*(1.0 -rect)*(Nd*(Nd-1.0))*vol;
@ -46,7 +46,7 @@ namespace Grid{
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu); U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
WilsonLoops<GaugeField>::RectStapleDouble(U2[mu],U[mu],mu); WilsonLoops<Gimpl>::RectStapleDouble(U2[mu],U[mu],mu);
} }
GaugeLinkField dSdU_mu(grid); GaugeLinkField dSdU_mu(grid);
@ -56,13 +56,11 @@ namespace Grid{
// Staple in direction mu // Staple in direction mu
WilsonLoops<GaugeField>::Staple(staple,Umu,mu); WilsonLoops<Gimpl>::Staple(staple,Umu,mu);
dSdU_mu = Ta(U[mu]*staple)*factor_p; dSdU_mu = Ta(U[mu]*staple)*factor_p;
// WilsonLoops<GaugeField>::RectStaple(staple,Umu,mu); WilsonLoops<Gimpl>::RectStaple(Umu,staple,U2,U,mu);
WilsonLoops<GaugeField>::RectStapleOptimised(staple,U2,U,mu);
dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r; dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r;
@ -78,31 +76,35 @@ namespace Grid{
// RBC c1 parameterisation is not really RBC but don't have good // RBC c1 parameterisation is not really RBC but don't have good
// reference and we are happy to change name if prior use of this plaq coeff // reference and we are happy to change name if prior use of this plaq coeff
// parameterisation is made known to us. // parameterisation is made known to us.
template<class GaugeField> template<class Gimpl>
class RBCGaugeAction : public PlaqPlusRectangleAction<GaugeField> { class RBCGaugeAction : public PlaqPlusRectangleAction<Gimpl> {
public: public:
RBCGaugeAction(RealD beta,RealD c1) : PlaqPlusRectangleAction<GaugeField>(beta*(1.0-8.0*c1), beta*c1) { INHERIT_GIMPL_TYPES(Gimpl);
RBCGaugeAction(RealD beta,RealD c1) : PlaqPlusRectangleAction<Gimpl>(beta*(1.0-8.0*c1), beta*c1) {
}; };
}; };
template<class GaugeField> template<class Gimpl>
class IwasakiGaugeAction : public RBCGaugeAction<GaugeField> { class IwasakiGaugeAction : public RBCGaugeAction<Gimpl> {
public: public:
IwasakiGaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-0.331) { INHERIT_GIMPL_TYPES(Gimpl);
IwasakiGaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-0.331) {
}; };
}; };
template<class GaugeField> template<class Gimpl>
class SymanzikGaugeAction : public RBCGaugeAction<GaugeField> { class SymanzikGaugeAction : public RBCGaugeAction<Gimpl> {
public: public:
SymanzikGaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-1.0/12.0) { INHERIT_GIMPL_TYPES(Gimpl);
SymanzikGaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-1.0/12.0) {
}; };
}; };
template<class GaugeField> template<class Gimpl>
class DBW2GaugeAction : public RBCGaugeAction<GaugeField> { class DBW2GaugeAction : public RBCGaugeAction<Gimpl> {
public: public:
DBW2GaugeAction(RealD beta) : RBCGaugeAction<GaugeField>(beta,-1.4067) { INHERIT_GIMPL_TYPES(Gimpl);
DBW2GaugeAction(RealD beta) : RBCGaugeAction<Gimpl>(beta,-1.4067) {
}; };
}; };

12
lib/qcd/action/gauge/WilsonGaugeAction.h
View File

@ -7,11 +7,13 @@ namespace Grid{
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Wilson Gauge Action .. should I template the Nc etc.. // Wilson Gauge Action .. should I template the Nc etc..
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
template<class GaugeField> template<class Gimpl>
class WilsonGaugeAction : public Action<GaugeField> { class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
public: public:
typedef LorentzScalar<GaugeField> GaugeLinkField; INHERIT_GIMPL_TYPES(Gimpl);
// typedef LorentzScalar<GaugeField> GaugeLinkField;
private: private:
RealD beta; RealD beta;
@ -21,7 +23,7 @@ namespace Grid{
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U); RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
RealD vol = U._grid->gSites(); RealD vol = U._grid->gSites();
RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5; RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
return action; return action;
@ -41,7 +43,7 @@ namespace Grid{
Umu = PeekIndex<LorentzIndex>(U,mu); Umu = PeekIndex<LorentzIndex>(U,mu);
// Staple in direction mu // Staple in direction mu
WilsonLoops<GaugeField>::Staple(dSdU_mu,U,mu); WilsonLoops<Gimpl>::Staple(dSdU_mu,U,mu);
dSdU_mu = Ta(Umu*dSdU_mu)*factor; dSdU_mu = Ta(Umu*dSdU_mu)*factor;
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu); PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
} }

20
lib/qcd/hmc/HMC.h
View File

@ -39,11 +39,12 @@ namespace Grid{
virtual void TrajectoryComplete (int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )=0; virtual void TrajectoryComplete (int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )=0;
}; };
template<class GaugeField> template<class Gimpl>
class PlaquetteLogger : public HmcObservable<GaugeField> { class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> {
private: private:
std::string Stem; std::string Stem;
public: public:
INHERIT_GIMPL_TYPES(Gimpl);
PlaquetteLogger(std::string cf) { PlaquetteLogger(std::string cf) {
Stem = cf; Stem = cf;
}; };
@ -52,11 +53,16 @@ namespace Grid{
{ {
std::string file; { std::ostringstream os; os << Stem <<"."<< traj; file = os.str(); } std::string file; { std::ostringstream os; os << Stem <<"."<< traj; file = os.str(); }
std::ofstream of(file); std::ofstream of(file);
RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
RealD rect = WilsonLoops<GaugeField>::avgRectangle(U); RealD peri_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(U);
of << plaq << " " << rect << std::endl; RealD peri_rect = WilsonLoops<PeriodicGimplR>::avgRectangle(U);
std::cout<< GridLogMessage<< "Plaquette for trajectory "<< traj << " is " << plaq <<std::endl;
std::cout<< GridLogMessage<< "Rectangle for trajectory "<< traj << " is " << rect <<std::endl; RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
RealD impl_rect = WilsonLoops<Gimpl>::avgRectangle(U);
of << traj<<" "<< impl_plaq << " " << impl_rect << " "<< peri_plaq<<" "<<peri_rect<<std::endl;
std::cout<< GridLogMessage<< "traj"<<" "<< "plaq " << " " << " rect " << " "<< "peri_plaq" <<" "<<"peri_rect"<<std::endl;
std::cout<< GridLogMessage<< traj<<" "<< impl_plaq << " " << impl_rect << " "<< peri_plaq<<" "<<peri_rect<<std::endl;
} }
}; };

28
lib/qcd/hmc/HmcRunner.h
View File

@ -5,12 +5,15 @@ namespace Grid{
namespace QCD{ namespace QCD{
class NerscHmcRunner { template<class Gimpl>
class NerscHmcRunnerTemplate {
public: public:
INHERIT_GIMPL_TYPES(Gimpl);
enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart }; enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
ActionSet<LatticeGaugeField> TheAction; ActionSet<GaugeField> TheAction;
GridCartesian * UGrid ; GridCartesian * UGrid ;
GridCartesian * FGrid ; GridCartesian * FGrid ;
@ -52,13 +55,13 @@ public:
} }
// Create integrator // Create integrator
typedef MinimumNorm2<LatticeGaugeField> IntegratorType;// change here to change the algorithm typedef MinimumNorm2<GaugeField> IntegratorType;// change here to change the algorithm
IntegratorParameters MDpar(20); IntegratorParameters MDpar(20);
IntegratorType MDynamics(UGrid,MDpar, TheAction); IntegratorType MDynamics(UGrid,MDpar, TheAction);
// Checkpoint strategy // Checkpoint strategy
NerscHmcCheckpointer<LatticeGaugeField> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1); NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
PlaquetteLogger<LatticeGaugeField> PlaqLog(std::string("plaq")); PlaquetteLogger<Gimpl> PlaqLog(std::string("plaq"));
HMCparameters HMCpar; HMCparameters HMCpar;
HMCpar.StartTrajectory = StartTraj; HMCpar.StartTrajectory = StartTraj;
@ -99,7 +102,7 @@ public:
Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG); Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
} }
HybridMonteCarlo<LatticeGaugeField,IntegratorType> HMC(HMCpar, MDynamics,sRNG,pRNG,U); HybridMonteCarlo<GaugeField,IntegratorType> HMC(HMCpar, MDynamics,sRNG,pRNG,U);
HMC.AddObservable(&Checkpoint); HMC.AddObservable(&Checkpoint);
HMC.AddObservable(&PlaqLog); HMC.AddObservable(&PlaqLog);
@ -109,5 +112,18 @@ public:
} }
}; };
typedef NerscHmcRunnerTemplate<PeriodicGimplR> NerscHmcRunner;
typedef NerscHmcRunnerTemplate<PeriodicGimplF> NerscHmcRunnerF;
typedef NerscHmcRunnerTemplate<PeriodicGimplD> NerscHmcRunnerD;
typedef NerscHmcRunnerTemplate<PeriodicGimplR> PeriodicNerscHmcRunner;
typedef NerscHmcRunnerTemplate<PeriodicGimplF> PeriodicNerscHmcRunnerF;
typedef NerscHmcRunnerTemplate<PeriodicGimplD> PeriodicNerscHmcRunnerD;
typedef NerscHmcRunnerTemplate<ConjugateGimplR> ConjugateNerscHmcRunner;
typedef NerscHmcRunnerTemplate<ConjugateGimplF> ConjugateNerscHmcRunnerF;
typedef NerscHmcRunnerTemplate<ConjugateGimplD> ConjugateNerscHmcRunnerD;
}} }}
#endif #endif

6
lib/qcd/hmc/NerscCheckpointer.h
View File

@ -10,13 +10,15 @@ namespace Grid{
namespace QCD{ namespace QCD{
template<class GaugeField> template<class Gimpl>
class NerscHmcCheckpointer : public HmcObservable<GaugeField> { class NerscHmcCheckpointer : public HmcObservable<typename Gimpl::GaugeField> {
private: private:
std::string configStem; std::string configStem;
std::string rngStem; std::string rngStem;
int SaveInterval; int SaveInterval;
public: public:
INHERIT_GIMPL_TYPES(Gimpl);
NerscHmcCheckpointer(std::string cf, std::string rn,int savemodulo) { NerscHmcCheckpointer(std::string cf, std::string rn,int savemodulo) {
configStem = cf; configStem = cf;
rngStem = rn; rngStem = rn;

97
lib/qcd/utils/CovariantCshift.h
View File

@ -1,25 +1,102 @@
#ifndef QCD_UTILS_COVARIANT_CSHIFT_H #ifndef QCD_UTILS_COVARIANT_CSHIFT_H
#define QCD_UTILS_COVARIANT_CSHIFT_H #define QCD_UTILS_COVARIANT_CSHIFT_H
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Low performance implementation of CovariantCshift API // Low performance implementation of CovariantCshift API
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// Make these members of an Impl class for BC's. // Make these members of an Impl class for BC's.
template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link,
namespace PeriodicBC {
template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link,
int mu, int mu,
const Lattice<covariant> &field) const Lattice<covariant> &field)
{ {
return Link*Cshift(field,mu,1);// moves towards negative mu return Link*Cshift(field,mu,1);// moves towards negative mu
}
template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field)
{
Lattice<covariant> tmp(field._grid);
tmp = adj(Link)*field;
return Cshift(tmp,mu,-1);// moves towards positive mu
}
} }
template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field) namespace GparityBC {
{
Lattice<covariant> tmp(field._grid); // Must give right answers across boundary
tmp = adj(Link)*field; // <----
return Cshift(tmp,mu,-1);// moves towards positive mu // --
// | |
// xxxxxxxxxxxxxxxxxxxx
// | |
//
// Stap= Cshift(GImpl::CovShiftForward(U[nu],nu,
// GImpl::CovShiftForward(U[nu],nu,
// GImpl::CovShiftBackward(U[mu],mu,
// GImpl::CovShiftBackward(U[nu],nu,
// GImpl::CovShiftIdentityBackward(U[nu],nu,-1))))) , mu, 1);
//
// U U^* U^* U^T U^adj = U (U U U^dag U^T )^*
// = U (U U U^dag)^* ( U^T )^*
//
// So covariant shift rule: conjugate inward shifted plane when crossing boundary applies.
//
// This conjugate should be applied to BOTH the link and the covariant field on backward shift
// boundary wrap.
//
// | |
// xxxxxxxxxxxxxxxxx
// | | <---- this link is conjugated, and the path leading into it. Segment crossing in and out is double conjugated.
// --
// ------->
template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field)
{
GridBase * grid = Link._grid;
int Lmu = grid->GlobalDimensions()[mu]-1;
conformable(field,Link);
Lattice<iScalar<vInteger> > coor(grid); LatticeCoordinate(coor,mu);
Lattice<covariant> field_bc = Cshift(field,mu,1);// moves towards negative mu;
field_bc = where(coor==Lmu,conjugate(field_bc),field_bc);
// std::cout<<"Gparity::CovCshiftForward mu="<<mu<<std::endl;
return Link*field_bc;
}
template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field)
{
GridBase * grid = field._grid;
int Lmu = grid->GlobalDimensions()[mu]-1;
conformable(field,Link);
Lattice<iScalar<vInteger> > coor(grid); LatticeCoordinate(coor,mu);
Lattice<covariant> tmp(grid);
tmp = adj(Link)*field;
tmp = where(coor==Lmu,conjugate(tmp),tmp);
// std::cout<<"Gparity::CovCshiftBackward mu="<<mu<<std::endl;
return Cshift(tmp,mu,-1);// moves towards positive mu
}
} }
}} }}
#endif #endif

173
lib/qcd/utils/WilsonLoops.h
View File

@ -4,18 +4,27 @@ namespace Grid {
namespace QCD { namespace QCD {
// Common wilson loop observables // Common wilson loop observables
template<class GaugeLorentz> template<class Gimpl>
class WilsonLoops { class WilsonLoops : public Gimpl {
public: public:
INHERIT_GIMPL_TYPES(Gimpl);
typedef LorentzScalar<GaugeLorentz> GaugeMat; typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
////////////////////////////////////////////////// //////////////////////////////////////////////////
// directed plaquette oriented in mu,nu plane // directed plaquette oriented in mu,nu plane
////////////////////////////////////////////////// //////////////////////////////////////////////////
static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu) static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
{ {
plaq=CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu])); // Annoyingly, must use either scope resolution to find dependent base class,
// or this-> ; there is no "this" in a static method. This forces explicit Gimpl scope
// resolution throughout the usage in this file, and rather defeats the purpose of deriving
// from Gimpl.
plaq= Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftBackward(U[nu],nu,
Gimpl::CovShiftForward (U[mu],mu,U[nu])));
} }
////////////////////////////////////////////////// //////////////////////////////////////////////////
// trace of directed plaquette oriented in mu,nu plane // trace of directed plaquette oriented in mu,nu plane
@ -98,10 +107,10 @@ public:
for(int d=0;d<Nd;d++){ for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d); U[d] = PeekIndex<LorentzIndex>(Umu,d);
} }
staple = zero; staple = zero;
GaugeMat tmp(grid); GaugeMat tmp(grid);
for(int nu=0;nu<Nd;nu++){ for(int nu=0;nu<Nd;nu++){
if(nu != mu) { if(nu != mu) {
@ -110,35 +119,24 @@ public:
// ^ // ^
// |__> nu // |__> nu
// __ __ // __
// | | // |
// __| = __| * // __|
// //
// staple += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1); staple+=Gimpl::ShiftStaple(
staple += Gimpl::CovShiftForward (U[nu],nu,
Cshift(CovShiftForward (U[nu],nu, Gimpl::CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,Cshift(adj(U[nu]),nu,-1))),mu,1); Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu);
// __
// Unu(x) Umu(x+nu) Unu^dag(x+mu) ; left mult by Umu^dag(x) to close ring // |
// |__
//
// __ __
// | |
// |__ = | * __
// //
// //
staple += staple+=Gimpl::ShiftStaple(
Cshift( Gimpl::CovShiftBackward(U[nu],nu,
CovShiftBackward(U[nu],nu, Gimpl::CovShiftBackward(U[mu],mu,U[nu])),mu);
CovShiftBackward(U[mu],mu,U[nu])),mu,1);
// tmp = CovShiftForward (U[mu],mu,U[nu]);
// staple+= CovShiftBackward(U[nu],nu,tmp);
// Unu^dag(x-nu) Umu(x-nu) Unu(x+mu-nu) ; left mult by Umu^dag(x) to close ring.
} }
} }
} }
@ -148,11 +146,11 @@ public:
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu) static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu)
{ {
rect = CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->| rect = Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) ; adj(Gimpl::CovShiftForward(U[nu],nu,Gimpl::CovShiftForward(U[mu],mu,U[mu]))) ;
rect = rect + rect = rect +
CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->|| Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[nu],nu,U[nu]))* // ->||
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) ; adj(Gimpl::CovShiftForward(U[nu],nu,Gimpl::CovShiftForward(U[nu],nu,U[mu]))) ;
} }
static void traceDirRectangle(LatticeComplex &rect, const std::vector<GaugeMat> &U, const int mu, const int nu) static void traceDirRectangle(LatticeComplex &rect, const std::vector<GaugeMat> &U, const int mu, const int nu)
{ {
@ -209,6 +207,12 @@ public:
static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){ static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){
U2 = U * Cshift(U,mu,1); U2 = U * Cshift(U,mu,1);
} }
////////////////////////////////////////////////////////////////////////////
// Hop by two optimisation strategy does not work nicely with Gparity. (could do,
// but need to track two deep where cross boundary and apply a conjugation).
// Must differentiate this in Gimpl, and use Gimpl::isPeriodicGaugeField to do so .
////////////////////////////////////////////////////////////////////////////
static void RectStapleOptimised(GaugeMat &Stap,std::vector<GaugeMat> &U2,std::vector<GaugeMat> &U,int mu){ static void RectStapleOptimised(GaugeMat &Stap,std::vector<GaugeMat> &U2,std::vector<GaugeMat> &U,int mu){
Stap = zero; Stap = zero;
@ -227,14 +231,14 @@ public:
tmp = adj(U2[mu])*tmp; tmp = adj(U2[mu])*tmp;
tmp = Cshift(tmp,mu,-2); tmp = Cshift(tmp,mu,-2);
Staple2x1 = CovShiftForward (U[nu],nu,tmp); Staple2x1 = Gimpl::CovShiftForward (U[nu],nu,tmp);
// Down staple // Down staple
// |___ ___| // |___ ___|
// //
tmp = adj(U2[mu])*U[nu]; tmp = adj(U2[mu])*U[nu];
Staple2x1+= CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2)); Staple2x1+= Gimpl::CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
// ___ ___ // ___ ___
@ -242,7 +246,7 @@ public:
// |___ ___| // |___ ___|
// //
Stap+= Cshift(CovShiftForward (U[mu],mu,Staple2x1),mu,1); Stap+= Cshift(Gimpl::CovShiftForward (U[mu],mu,Staple2x1),mu,1);
// ___ ___ // ___ ___
// |___ | // |___ |
@ -259,7 +263,7 @@ public:
// | | // | |
tmp = Cshift(adj(U2[nu]),nu,-2); tmp = Cshift(adj(U2[nu]),nu,-2);
tmp = CovShiftBackward(U[mu],mu,tmp); tmp = Gimpl::CovShiftBackward(U[mu],mu,tmp);
tmp = U2[nu]*Cshift(tmp,nu,2); tmp = U2[nu]*Cshift(tmp,nu,2);
Stap+= Cshift(tmp, mu, 1); Stap+= Cshift(tmp, mu, 1);
@ -268,7 +272,7 @@ public:
// | | // | |
// -- // --
tmp = CovShiftBackward(U[mu],mu,U2[nu]); tmp = Gimpl::CovShiftBackward(U[mu],mu,U2[nu]);
tmp = adj(U2[nu])*tmp; tmp = adj(U2[nu])*tmp;
tmp = Cshift(tmp,nu,-2); tmp = Cshift(tmp,nu,-2);
Stap+=Cshift(tmp, mu, 1); Stap+=Cshift(tmp, mu, 1);
@ -276,23 +280,20 @@ public:
} }
static void RectStaple(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
#if 1 static void RectStaple(GaugeMat &Stap,const GaugeLorentz & Umu,int mu)
{
RectStapleUnoptimised(Stap,Umu,mu); RectStapleUnoptimised(Stap,Umu,mu);
#else }
GridBase *grid = Umu._grid; static void RectStaple(const GaugeLorentz & Umu,GaugeMat &Stap,
std::vector<GaugeMat> &U2,
std::vector<GaugeMat> U (Nd,grid); std::vector<GaugeMat> &U, int mu)
std::vector<GaugeMat> U2(Nd,grid); {
if ( Gimpl::isPeriodicGaugeField() ){
for(int mu=0;mu<Nd;mu++){ RectStapleOptimised(Stap,U2,U,mu);
U[mu] = PeekIndex<LorentzIndex>(Umu,mu); } else {
RectStapleDouble(U2[mu],U[mu],mu); RectStapleUnoptimised(Stap,Umu,mu);
} }
RectStapleOptimised(Stap,U2,U,mu);
#endif
} }
static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){ static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
@ -310,46 +311,51 @@ public:
// __ ___ // __ ___
// | __ | // | __ |
// //
Stap+= Cshift(CovShiftForward (U[mu],mu, Stap+= Gimpl::ShiftStaple(
CovShiftForward (U[nu],nu, Gimpl::CovShiftForward (U[mu],mu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[mu],mu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1); Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftIdentityBackward(U[nu],nu))))) , mu);
// __ // __
// |__ __ | // |__ __ |
Stap+= Cshift(CovShiftForward (U[mu],mu, Stap+= Gimpl::ShiftStaple(
CovShiftBackward(U[nu],nu, Gimpl::CovShiftForward (U[mu],mu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1); Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftBackward(U[mu],mu, U[nu])))) , mu);
// __ // __
// |__ __ | // |__ __ |
Stap+= Cshift(CovShiftBackward(U[nu],nu, Stap+= Gimpl::ShiftStaple(
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[mu],mu,
CovShiftForward(U[nu],nu,U[mu])))) , mu, 1); Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftForward(U[nu],nu,U[mu])))) , mu);
// __ ___ // __ ___
// |__ | // |__ |
Stap+= Cshift(CovShiftForward (U[nu],nu, Stap+= Gimpl::ShiftStaple(
CovShiftBackward(U[mu],mu, Gimpl::CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[mu],mu,
CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1); Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftBackward(U[nu],nu,U[mu])))) , mu);
// -- // --
// | | // | |
// //
// | | // | |
Stap+= Cshift(CovShiftForward(U[nu],nu, Stap+= Gimpl::ShiftStaple(
CovShiftForward(U[nu],nu, Gimpl::CovShiftForward(U[nu],nu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftForward(U[nu],nu,
CovShiftBackward(U[nu],nu, Gimpl::CovShiftBackward(U[mu],mu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1); Gimpl::CovShiftBackward(U[nu],nu,
Gimpl::CovShiftIdentityBackward(U[nu],nu))))) , mu);
// | | // | |
@ -357,10 +363,11 @@ public:
// | | // | |
// -- // --
Stap+= Cshift(CovShiftBackward(U[nu],nu, Stap+= Gimpl::ShiftStaple(
CovShiftBackward(U[nu],nu, Gimpl::CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(U[nu],nu,
CovShiftForward (U[nu],nu,U[nu])))) , mu, 1); Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftForward (U[nu],nu,U[nu])))) , mu);
}} }}
} }
@ -368,10 +375,10 @@ public:
}; };
typedef WilsonLoops<LatticeGaugeField> ColourWilsonLoops; typedef WilsonLoops<PeriodicGimplR> ColourWilsonLoops;
typedef WilsonLoops<LatticeGaugeField> U1WilsonLoops; typedef WilsonLoops<PeriodicGimplR> U1WilsonLoops;
typedef WilsonLoops<LatticeGaugeField> SU2WilsonLoops; typedef WilsonLoops<PeriodicGimplR> SU2WilsonLoops;
typedef WilsonLoops<LatticeGaugeField> SU3WilsonLoops; typedef WilsonLoops<PeriodicGimplR> SU3WilsonLoops;
}} }}

14
tests/Make.inc
View File

@ -1,5 +1,5 @@
bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gp_rect_force Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd
Test_GaugeAction_SOURCES=Test_GaugeAction.cc Test_GaugeAction_SOURCES=Test_GaugeAction.cc
@ -102,6 +102,10 @@ Test_gamma_SOURCES=Test_gamma.cc
Test_gamma_LDADD=-lGrid Test_gamma_LDADD=-lGrid
Test_gp_rect_force_SOURCES=Test_gp_rect_force.cc
Test_gp_rect_force_LDADD=-lGrid
Test_gparity_SOURCES=Test_gparity.cc Test_gparity_SOURCES=Test_gparity.cc
Test_gparity_LDADD=-lGrid Test_gparity_LDADD=-lGrid
@ -130,6 +134,14 @@ Test_hmc_EOWilsonRatio_SOURCES=Test_hmc_EOWilsonRatio.cc
Test_hmc_EOWilsonRatio_LDADD=-lGrid Test_hmc_EOWilsonRatio_LDADD=-lGrid
Test_hmc_GparityIwasakiGauge_SOURCES=Test_hmc_GparityIwasakiGauge.cc
Test_hmc_GparityIwasakiGauge_LDADD=-lGrid
Test_hmc_GparityWilsonGauge_SOURCES=Test_hmc_GparityWilsonGauge.cc
Test_hmc_GparityWilsonGauge_LDADD=-lGrid
Test_hmc_IwasakiGauge_SOURCES=Test_hmc_IwasakiGauge.cc Test_hmc_IwasakiGauge_SOURCES=Test_hmc_IwasakiGauge.cc
Test_hmc_IwasakiGauge_LDADD=-lGrid Test_hmc_IwasakiGauge_LDADD=-lGrid

2
tests/Test_GaugeAction.cc
View File

@ -87,7 +87,7 @@ int main (int argc, char ** argv)
Plaq = zero; Plaq = zero;
for(int mu=1;mu<Nd;mu++){ for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){ for(int nu=0;nu<mu;nu++){
Plaq = Plaq + trace(CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]))); Plaq = Plaq + trace(PeriodicBC::CovShiftForward(U[mu],mu,U[nu])*adj(PeriodicBC::CovShiftForward(U[nu],nu,U[mu])));
} }
} }

22
tests/Test_RectPlaq.cc
View File

@ -34,11 +34,11 @@ void RectPlaq(const std::vector<LatticeColourMatrix> &U, LatticeComplex &RectPla
for(int mu=1;mu<Nd;mu++){ for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){ for(int nu=0;nu<mu;nu++){
RectPlaqValue = RectPlaqValue + trace( RectPlaqValue = RectPlaqValue + trace(
CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->| PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) ); adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[mu],mu,U[mu]))) );
RectPlaqValue = RectPlaqValue + trace( RectPlaqValue = RectPlaqValue + trace(
CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->|| PeriodicBC::CovShiftForward(U[mu],mu,PeriodicBC::CovShiftForward(U[nu],nu,U[nu]))* // ->||
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) ); adj(PeriodicBC::CovShiftForward(U[nu],nu,PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) );
} }
} }
} }
@ -87,7 +87,7 @@ int main (int argc, char ** argv)
Plaq = zero; Plaq = zero;
for(int mu=1;mu<Nd;mu++){ for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){ for(int nu=0;nu<mu;nu++){
Plaq = Plaq + trace(CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]))); Plaq = Plaq + trace(PeriodicBC::CovShiftForward(U[mu],mu,U[nu])*adj(PeriodicBC::CovShiftForward(U[nu],nu,U[mu])));
} }
} }
@ -134,7 +134,7 @@ int main (int argc, char ** argv)
/* /*
(x) ---> ---> : U(x,mu)*U(x+mu, mu) (x) ---> ---> : U(x,mu)*U(x+mu, mu)
*/ */
left_2= CovShiftForward(U[mu],mu,U[mu]); left_2= PeriodicBC::CovShiftForward(U[mu],mu,U[mu]);
/* /*
upper_l = <---- <--- upper_l = <---- <---
^ ^
@ -144,7 +144,7 @@ int main (int argc, char ** argv)
*/ */
tmp=Cshift(U[nu],mu,2); tmp=Cshift(U[nu],mu,2);
upper_l= CovShiftForward(tmp,nu,adj(left_2)); // i.e. upper_l upper_l= PeriodicBC::CovShiftForward(tmp,nu,adj(left_2)); // i.e. upper_l
/* /*
upper_staple= <---- <--- ^ upper_staple= <---- <--- ^
| | | |
@ -170,7 +170,7 @@ int main (int argc, char ** argv)
| | | |
<-----<--- V <-----<--- V
*/ */
tmp=CovShiftBackward(U[mu],nu,down_staple); tmp=PeriodicBC::CovShiftBackward(U[mu],nu,down_staple);
ds_U+=Cshift(tmp,mu,-1); ds_U+=Cshift(tmp,mu,-1);
/* /*
ds_U+= <----<---- ^ ds_U+= <----<---- ^
@ -193,7 +193,7 @@ int main (int argc, char ** argv)
| |
(x) (x)
*/ */
LatticeColourMatrix up2= CovShiftForward(U[nu],nu,U[nu]); LatticeColourMatrix up2= PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
/* /*
<----^ <----^
| |
@ -202,7 +202,7 @@ int main (int argc, char ** argv)
(x) (x)
*/ */
// Unu(x+mu)Unu(x+mu+nu) UmuDag(x+nu+nu) lives at X // Unu(x+mu)Unu(x+mu+nu) UmuDag(x+nu+nu) lives at X
upper_l= CovShiftForward(Cshift(up2,mu,1),nu,Cshift(adj(U[mu]),nu,1)); upper_l= PeriodicBC::CovShiftForward(Cshift(up2,mu,1),nu,Cshift(adj(U[mu]),nu,1));
/* /*
|<----^ |<----^
upper_staple = V | upper_staple = V |
@ -218,7 +218,7 @@ int main (int argc, char ** argv)
downer_l= | downer_l= |
(x)<----V (x)<----V
*/ */
upper_l= adj(CovShiftForward(U[mu],mu,up2)); //downer_l upper_l= adj(PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
/* /*
^ | ^ |
down_staple = | V down_staple = | V

98
tests/Test_gp_rect_force.cc Normal file
View File

@ -0,0 +1,98 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedRandomDevice();
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
double beta = 1.0;
double c1 = 0.331;
//GparityPlaqPlusRectangleActionR Action(beta,c1);
ConjugateWilsonGaugeActionR Action(beta);
//WilsonGaugeActionR Action(beta);
ComplexD S = Action.S(U);
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(&Grid);
Action.deriv(U,UdSdU);
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
LatticeColourMatrix mommu(&Grid);
LatticeColourMatrix forcemu(&Grid);
LatticeGaugeField mom(&Grid);
LatticeGaugeField Uprime(&Grid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
// fourth order exponential approx
parallel_for(auto i=mom.begin();i<mom.end();i++){ // exp(pmu dt) * Umu
Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
}
}
ComplexD Sprime = Action.S(Uprime);
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(&Grid); dS = zero;
for(int mu=0;mu<Nd;mu++){
auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update gauge action density
// U = exp(p dt) U
// dU/dt = p U
// so dSdt = trace( dUdt dSdU) = trace( p UdSdUmu )
dS = dS - trace(mommu*UdSdUmu)*dt*2.0;
}
Complex dSpred = sum(dS);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "pred dS "<< dSpred <<std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

4
tests/Test_hmc_EODWFRatio_Gparity.cc
View File

@ -8,7 +8,7 @@ namespace Grid {
namespace QCD { namespace QCD {
class HmcRunner : public NerscHmcRunner { class HmcRunner : public ConjugateNerscHmcRunner {
public: public:
void BuildTheAction (int argc, char **argv) void BuildTheAction (int argc, char **argv)
@ -31,7 +31,7 @@ public:
LatticeGaugeField U(UGrid); LatticeGaugeField U(UGrid);
// Gauge action // Gauge action
WilsonGaugeActionR Waction(5.6); ConjugateWilsonGaugeActionR Waction(5.6);
Real mass=0.04; Real mass=0.04;
Real pv =1.0; Real pv =1.0;

53
tests/Test_hmc_GparityIwasakiGauge.cc Normal file
View File

@ -0,0 +1,53 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public ConjugateNerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
ConjugateIwasakiGaugeActionR Gaction(2.6);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Gaction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);
}

53
tests/Test_hmc_GparityWilsonGauge.cc Normal file
View File

@ -0,0 +1,53 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public ConjugateNerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
ConjugateWilsonGaugeActionR Waction(5.6);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Waction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);
}