portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-03 21:44:33 +00:00
Code Releases Activity

Cleaned up HMC output. Tested smeared HMCs for single precision (OK)

Browse Source
This commit is contained in:
Guido Cossu
2016-07-05 12:03:54 +01:00
parent fdfbf11c6d
commit 3e80947c2b
6 changed files with 838 additions and 787 deletions
Download Patch File Download Diff File Expand all files Collapse all files

101
lib/Log.cc
View File

@@ -1,51 +1,51 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/Log.cc Source file: ./lib/Log.cc
Copyright (C) 2015 Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk> Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
#include <Grid.h> #include <Grid.h>
namespace Grid { namespace Grid {
GridStopWatch Logger::StopWatch;
std::ostream Logger::devnull(0);
Colours GridLogColours (0);
GridLogger GridLogError (1,"Error",GridLogColours, "RED");
GridLogger GridLogWarning (1,"Warning",GridLogColours, "YELLOW");
GridLogger GridLogMessage (1,"Message",GridLogColours, "NORMAL");
GridLogger GridLogDebug (1,"Debug",GridLogColours, "PURPLE");
GridLogger GridLogPerformance(1,"Performance",GridLogColours, "GREEN");
GridLogger GridLogIterative (1,"Iterative",GridLogColours, "BLUE");
GridLogger GridLogIntegrator (1,"Integrator",GridLogColours, "BLUE");
void GridLogConfigure(std::vector<std::string> &logstreams) GridStopWatch Logger::StopWatch;
{ std::ostream Logger::devnull(0);
Colours GridLogColours(0);
GridLogger GridLogError(1, "Error", GridLogColours, "RED");
GridLogger GridLogWarning(1, "Warning", GridLogColours, "YELLOW");
GridLogger GridLogMessage(1, "Message", GridLogColours, "NORMAL");
GridLogger GridLogDebug(1, "Debug", GridLogColours, "PURPLE");
GridLogger GridLogPerformance(1, "Performance", GridLogColours, "GREEN");
GridLogger GridLogIterative(1, "Iterative", GridLogColours, "BLUE");
GridLogger GridLogIntegrator(1, "Integrator", GridLogColours, "BLUE");
void GridLogConfigure(std::vector<std::string> &logstreams) {
GridLogError.Active(0); GridLogError.Active(0);
GridLogWarning.Active(0); GridLogWarning.Active(0);
GridLogMessage.Active(0); GridLogMessage.Active(0);
@@ -55,43 +55,38 @@ void GridLogConfigure(std::vector<std::string> &logstreams)
GridLogIntegrator.Active(0); GridLogIntegrator.Active(0);
GridLogColours.Active(0); GridLogColours.Active(0);
for(int i=0;i<logstreams.size();i++){ for (int i = 0; i < logstreams.size(); i++) {
if ( logstreams[i]== std::string("Error") ) GridLogError.Active(1); if (logstreams[i] == std::string("Error")) GridLogError.Active(1);
if ( logstreams[i]== std::string("Warning") ) GridLogWarning.Active(1); if (logstreams[i] == std::string("Warning")) GridLogWarning.Active(1);
if ( logstreams[i]== std::string("Message") ) GridLogMessage.Active(1); if (logstreams[i] == std::string("Message")) GridLogMessage.Active(1);
if ( logstreams[i]== std::string("Iterative") ) GridLogIterative.Active(1); if (logstreams[i] == std::string("Iterative")) GridLogIterative.Active(1);
if ( logstreams[i]== std::string("Debug") ) GridLogDebug.Active(1); if (logstreams[i] == std::string("Debug")) GridLogDebug.Active(1);
if ( logstreams[i]== std::string("Performance") ) GridLogPerformance.Active(1); if (logstreams[i] == std::string("Performance"))
if ( logstreams[i]== std::string("Integrator" ) ) GridLogIntegrator.Active(1); GridLogPerformance.Active(1);
if ( logstreams[i]== std::string("Colours" ) ) GridLogColours.Active(1); if (logstreams[i] == std::string("Integrator")) GridLogIntegrator.Active(1);
if (logstreams[i] == std::string("Colours")) GridLogColours.Active(1);
} }
} }
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
// Verbose limiter on MPI tasks // Verbose limiter on MPI tasks
//////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////
void Grid_quiesce_nodes(void) void Grid_quiesce_nodes(void) {
{ int me = 0;
int me=0;
#ifdef GRID_COMMS_MPI #ifdef GRID_COMMS_MPI
MPI_Comm_rank(MPI_COMM_WORLD,&me); MPI_Comm_rank(MPI_COMM_WORLD, &me);
#endif #endif
#ifdef GRID_COMMS_SHMEM #ifdef GRID_COMMS_SHMEM
me = shmem_my_pe(); me = shmem_my_pe();
#endif #endif
if ( me ) { if (me) {
std::cout.setstate(std::ios::badbit); std::cout.setstate(std::ios::badbit);
} }
} }
void Grid_unquiesce_nodes(void) void Grid_unquiesce_nodes(void) {
{
#ifdef GRID_COMMS_MPI #ifdef GRID_COMMS_MPI
std::cout.clear(); std::cout.clear();
#endif #endif
} }
} }

7
lib/Log.h
View File

@@ -6,9 +6,9 @@
Copyright (C) 2015 Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk> Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@@ -49,7 +49,6 @@ protected:
public: public:
std::map<std::string, std::string> colour; std::map<std::string, std::string> colour;
Colours(bool activate=false){ Colours(bool activate=false){
Active(activate); Active(activate);
}; };

331
lib/qcd/action/gauge/GaugeImpl.h
View File

@@ -1,189 +1,188 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/gauge/GaugeImpl.h Source file: ./lib/qcd/action/gauge/GaugeImpl.h
Copyright (C) 2015 Copyright (C) 2015
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
#ifndef GRID_QCD_GAUGE_IMPL_H /* END LEGAL */
#define GRID_QCD_GAUGE_IMPL_H #ifndef GRID_QCD_GAUGE_IMPL_H
#define GRID_QCD_GAUGE_IMPL_H
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////// // Implementation dependent gauge types
// Implementation dependent gauge types ////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
template<class Gimpl> class WilsonLoops; template <class Gimpl> class WilsonLoops;
#define INHERIT_GIMPL_TYPES(GImpl) \ #define INHERIT_GIMPL_TYPES(GImpl) \
typedef typename GImpl::Simd Simd;\ typedef typename GImpl::Simd Simd; \
typedef typename GImpl::GaugeLinkField GaugeLinkField;\ typedef typename GImpl::GaugeLinkField GaugeLinkField; \
typedef typename GImpl::GaugeField GaugeField;\ typedef typename GImpl::GaugeField GaugeField; \
typedef typename GImpl::SiteGaugeField SiteGaugeField;\ typedef typename GImpl::SiteGaugeField SiteGaugeField; \
typedef typename GImpl::SiteGaugeLink SiteGaugeLink; typedef typename GImpl::SiteGaugeLink SiteGaugeLink;
// //
template<class S,int Nrepresentation=Nc> template <class S, int Nrepresentation = Nc> class GaugeImplTypes {
class GaugeImplTypes { public:
public: typedef S Simd;
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;
// Move this elsewhere? template <typename vtype>
static inline void AddGaugeLink(GaugeField& U, GaugeLinkField& W, int mu){ // U[mu] += W 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;
// Move this elsewhere?
static inline void AddGaugeLink(GaugeField &U, GaugeLinkField &W,
int mu) { // U[mu] += W
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(auto ss=0;ss<U._grid->oSites();ss++){ for (auto ss = 0; ss < U._grid->oSites(); ss++) {
U._odata[ss]._internal[mu] = U._odata[ss]._internal[mu] + W._odata[ss]._internal; U._odata[ss]._internal[mu] =
} U._odata[ss]._internal[mu] + W._odata[ss]._internal;
} }
};
// Composition with smeared link, bc's etc.. probably need multiple inheritance
// Variable precision "S" and variable Nc
template<class GimplTypes>
class PeriodicGaugeImpl : public GimplTypes {
public:
INHERIT_GIMPL_TYPES(GimplTypes);
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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 GimplTypes>
class ConjugateGaugeImpl : public GimplTypes {
public:
INHERIT_GIMPL_TYPES(GimplTypes);
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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 ConjugateBC::CovShiftForward(Link,mu,field);
}
template<class covariant> static
Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,const Lattice<covariant> &field) {
return ConjugateBC::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 GaugeImplTypes<vComplex,Nc> GimplTypesR;
typedef GaugeImplTypes<vComplexF,Nc> GimplTypesF;
typedef GaugeImplTypes<vComplexD,Nc> GimplTypesD;
typedef PeriodicGaugeImpl<GimplTypesR> PeriodicGimplR; // Real.. whichever prec
typedef PeriodicGaugeImpl<GimplTypesF> PeriodicGimplF; // Float
typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; // Double
typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever prec
typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
} }
};
// Composition with smeared link, bc's etc.. probably need multiple inheritance
// Variable precision "S" and variable Nc
template <class GimplTypes> class PeriodicGaugeImpl : public GimplTypes {
public:
INHERIT_GIMPL_TYPES(GimplTypes);
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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 GimplTypes> class ConjugateGaugeImpl : public GimplTypes {
public:
INHERIT_GIMPL_TYPES(GimplTypes);
////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 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 ConjugateBC::CovShiftForward(Link, mu, field);
}
template <class covariant>
static Lattice<covariant> CovShiftBackward(const GaugeLinkField &Link, int mu,
const Lattice<covariant> &field) {
return ConjugateBC::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 GaugeImplTypes<vComplex, Nc> GimplTypesR;
typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
typedef PeriodicGaugeImpl<GimplTypesR> PeriodicGimplR; // Real.. whichever prec
typedef PeriodicGaugeImpl<GimplTypesF> PeriodicGimplF; // Float
typedef PeriodicGaugeImpl<GimplTypesD> PeriodicGimplD; // Double
typedef ConjugateGaugeImpl<GimplTypesR>
ConjugateGimplR; // Real.. whichever prec
typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
}
} }
#endif #endif

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

@@ -1,33 +1,34 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/HMC.h Source file: ./lib/qcd/hmc/HMC.h
Copyright (C) 2015 Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk> Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: neo <cossu@post.kek.jp> Author: neo <cossu@post.kek.jp>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. (at your option) any later version.
This program is distributed in the hope that it will be useful, This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. GNU General Public License for more details.
You should have received a copy of the GNU General Public License along 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., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
*************************************************************************************/ directory
/* END LEGAL */ *************************************************************************************/
/* END LEGAL */
//-------------------------------------------------------------------- //--------------------------------------------------------------------
/*! @file HMC.h /*! @file HMC.h
* @brief Classes for Hybrid Monte Carlo update * @brief Classes for Hybrid Monte Carlo update
@@ -41,172 +42,195 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <string> #include <string>
namespace Grid {
namespace QCD {
namespace Grid{ struct HMCparameters {
namespace QCD{ Integer StartTrajectory;
Integer Trajectories; /* @brief Number of sweeps in this run */
bool MetropolisTest;
Integer NoMetropolisUntil;
struct HMCparameters{ HMCparameters() {
////////////////////////////// Default values
MetropolisTest = true;
NoMetropolisUntil = 10;
StartTrajectory = 0;
Trajectories = 200;
/////////////////////////////////
}
Integer StartTrajectory; void print() const {
Integer Trajectories; /* @brief Number of sweeps in this run */ std::cout << GridLogMessage << "[HMC parameter] Trajectories : " << Trajectories << "\n";
bool MetropolisTest; std::cout << GridLogMessage << "[HMC parameter] Start trajectory : " << StartTrajectory << "\n";
Integer NoMetropolisUntil; std::cout << GridLogMessage << "[HMC parameter] Metropolis test (on/off): " << MetropolisTest << "\n";
std::cout << GridLogMessage << "[HMC parameter] Thermalization trajs : " << NoMetropolisUntil << "\n";
}
};
HMCparameters(){ template <class GaugeField>
////////////////////////////// Default values class HmcObservable {
MetropolisTest = true; public:
NoMetropolisUntil = 10; virtual void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
StartTrajectory = 0; GridParallelRNG &pRNG) = 0;
Trajectories = 200; };
/////////////////////////////////
}
};
template<class GaugeField> template <class Gimpl>
class HmcObservable { class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> {
public: private:
virtual void TrajectoryComplete (int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )=0; std::string Stem;
};
template<class Gimpl> public:
class PlaquetteLogger : public HmcObservable<typename Gimpl::GaugeField> { INHERIT_GIMPL_TYPES(Gimpl);
private: PlaquetteLogger(std::string cf) { Stem = cf; };
std::string Stem;
public:
INHERIT_GIMPL_TYPES(Gimpl);
PlaquetteLogger(std::string cf) {
Stem = cf;
};
void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG ) void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
{ GridParallelRNG &pRNG) {
std::string file; { std::ostringstream os; os << Stem <<"."<< traj; file = os.str(); } std::string file;
std::ofstream of(file); {
std::ostringstream os;
os << Stem << "." << traj;
file = os.str();
}
std::ofstream of(file);
RealD peri_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(U); RealD peri_plaq = WilsonLoops<PeriodicGimplR>::avgPlaquette(U);
RealD peri_rect = WilsonLoops<PeriodicGimplR>::avgRectangle(U); RealD peri_rect = WilsonLoops<PeriodicGimplR>::avgRectangle(U);
RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(U); RealD impl_plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
RealD impl_rect = WilsonLoops<Gimpl>::avgRectangle(U); RealD impl_rect = WilsonLoops<Gimpl>::avgRectangle(U);
of << traj<<" "<< impl_plaq << " " << impl_rect << " "<< peri_plaq<<" "<<peri_rect<<std::endl; of << traj << " " << impl_plaq << " " << impl_rect << " " << peri_plaq
std::cout<< GridLogMessage<< "traj"<<" "<< "plaq " << " " << " rect " << " "<< "peri_plaq" <<" "<<"peri_rect"<<std::endl; << " " << peri_rect << std::endl;
std::cout<< GridLogMessage<< 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;
}
};
// template <class GaugeField, class Integrator, class Smearer, class Boundary> // template <class GaugeField, class Integrator, class Smearer, class
template <class GaugeField, class IntegratorType> // Boundary>
class HybridMonteCarlo { template <class GaugeField, class IntegratorType>
private: class HybridMonteCarlo {
private:
const HMCparameters Params;
const HMCparameters Params; GridSerialRNG &sRNG; // Fixme: need a RNG management strategy.
GridParallelRNG &pRNG; // Fixme: need a RNG management strategy.
GridSerialRNG &sRNG; // Fixme: need a RNG management strategy. GaugeField &Ucur;
GridParallelRNG &pRNG; // Fixme: need a RNG management strategy.
GaugeField & Ucur;
IntegratorType &TheIntegrator; IntegratorType &TheIntegrator;
std::vector<HmcObservable<GaugeField> *> Observables; std::vector<HmcObservable<GaugeField> *> Observables;
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
// Metropolis step // Metropolis step
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
bool metropolis_test(const RealD DeltaH){ bool metropolis_test(const RealD DeltaH) {
RealD rn_test;
RealD rn_test; RealD prob = std::exp(-DeltaH);
RealD prob = std::exp(-DeltaH); random(sRNG, rn_test);
random(sRNG,rn_test); std::cout << GridLogMessage
<< "--------------------------------------------------\n";
std::cout<<GridLogMessage<< "--------------------------------------------\n"; std::cout << GridLogMessage << "exp(-dH) = " << prob
std::cout<<GridLogMessage<< "dH = "<<DeltaH << " Random = "<< rn_test <<"\n"; << " Random = " << rn_test << "\n";
std::cout<<GridLogMessage<< "Acc. Probability = " << ((prob<1.0)? prob: 1.0)<< " "; std::cout << GridLogMessage
<< "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
if((prob >1.0) || (rn_test <= prob)){ // accepted
std::cout<<GridLogMessage <<"-- ACCEPTED\n";
return true;
} else { // rejected
std::cout<<GridLogMessage <<"-- REJECTED\n";
return false;
}
if ((prob > 1.0) || (rn_test <= prob)) { // accepted
std::cout << GridLogMessage << "Metropolis_test -- ACCEPTED\n";
std::cout << GridLogMessage
<< "--------------------------------------------------\n";
return true;
} else { // rejected
std::cout << GridLogMessage << "Metropolis_test -- REJECTED\n";
std::cout << GridLogMessage
<< "--------------------------------------------------\n";
return false;
}
}
/////////////////////////////////////////////////////////
// Evolution
/////////////////////////////////////////////////////////
RealD evolve_step(GaugeField &U) {
TheIntegrator.refresh(U, pRNG); // set U and initialize P and phi's
RealD H0 = TheIntegrator.S(U); // initial state action
std::streamsize current_precision = std::cout.precision();
std::cout.precision(17);
std::cout << GridLogMessage << "Total H before trajectory = " << H0 << "\n";
std::cout.precision(current_precision);
TheIntegrator.integrate(U);
RealD H1 = TheIntegrator.S(U); // updated state action
std::cout.precision(17);
std::cout << GridLogMessage << "Total H after trajectory = " << H1
<< " dH = " << H1 - H0 << "\n";
std::cout.precision(current_precision);
return (H1 - H0);
}
public:
/////////////////////////////////////////
// Constructor
/////////////////////////////////////////
HybridMonteCarlo(HMCparameters Pams, IntegratorType &_Int,
GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U)
: Params(Pams), TheIntegrator(_Int), sRNG(_sRNG), pRNG(_pRNG), Ucur(_U) {}
~HybridMonteCarlo(){};
void AddObservable(HmcObservable<GaugeField> *obs) {
Observables.push_back(obs);
}
void evolve(void) {
Real DeltaH;
GaugeField Ucopy(Ucur._grid);
Params.print();
// Actual updates (evolve a copy Ucopy then copy back eventually)
for (int traj = Params.StartTrajectory;
traj < Params.Trajectories + Params.StartTrajectory; ++traj) {
std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
Ucopy = Ucur;
DeltaH = evolve_step(Ucopy);
bool accept = true;
if (traj >= Params.NoMetropolisUntil) {
accept = metropolis_test(DeltaH);
} }
///////////////////////////////////////////////////////// if (accept) {
// Evolution Ucur = Ucopy;
/////////////////////////////////////////////////////////
RealD evolve_step(GaugeField& U){
TheIntegrator.refresh(U,pRNG); // set U and initialize P and phi's
RealD H0 = TheIntegrator.S(U); // initial state action
std::cout<<GridLogMessage<<"Total H before = "<< H0 << "\n";
TheIntegrator.integrate(U);
RealD H1 = TheIntegrator.S(U); // updated state action
std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
return (H1-H0);
}
public:
/////////////////////////////////////////
// Constructor
/////////////////////////////////////////
HybridMonteCarlo(HMCparameters Pms, IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U) :
Params(Pms),
TheIntegrator(_Int),
sRNG(_sRNG),
pRNG(_pRNG),
Ucur(_U)
{
}
~HybridMonteCarlo(){};
void AddObservable(HmcObservable<GaugeField> *obs) {
Observables.push_back(obs);
} }
void evolve(void){ for (int obs = 0; obs < Observables.size(); obs++) {
Observables[obs]->TrajectoryComplete(traj + 1, Ucur, sRNG, pRNG);
Real DeltaH;
GaugeField Ucopy(Ucur._grid);
// Actual updates (evolve a copy Ucopy then copy back eventually)
for(int traj=Params.StartTrajectory; traj < Params.Trajectories+Params.StartTrajectory; ++traj){
std::cout<<GridLogMessage << "-- # Trajectory = "<< traj << "\n";
Ucopy = Ucur;
DeltaH = evolve_step(Ucopy);
bool accept = true;
if ( traj > Params.NoMetropolisUntil) {
accept = metropolis_test(DeltaH);
}
if ( accept ) {
Ucur = Ucopy;
}
for(int obs = 0;obs<Observables.size();obs++){
Observables[obs]->TrajectoryComplete (traj+1,Ucur,sRNG,pRNG);
}
}
} }
}; }
}
}// QCD };
}// Grid
} // QCD
} // Grid
#endif #endif

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

@@ -81,6 +81,14 @@ public:
NumTraj = ivec[0]; NumTraj = ivec[0];
} }
int NumThermalizations = 10;
if( GridCmdOptionExists(argv,argv+argc,"--Thermalizations") ){
arg= GridCmdOptionPayload(argv,argv+argc,"--Thermalizations");
std::vector<int> ivec(0);
GridCmdOptionIntVector(arg,ivec);
NumThermalizations = ivec[0];
}
GridSerialRNG sRNG; GridSerialRNG sRNG;
GridParallelRNG pRNG(UGrid); GridParallelRNG pRNG(UGrid);
@@ -110,33 +118,30 @@ public:
PlaquetteLogger<Gimpl> PlaqLog(std::string("plaq")); PlaquetteLogger<Gimpl> PlaqLog(std::string("plaq"));
HMCparameters HMCpar; HMCparameters HMCpar;
HMCpar.StartTrajectory = StartTraj; HMCpar.StartTrajectory = StartTraj;
HMCpar.Trajectories = NumTraj; HMCpar.Trajectories = NumTraj;
HMCpar.NoMetropolisUntil = NumThermalizations;
if ( StartType == HotStart ) { if ( StartType == HotStart ) {
// Hot start // Hot start
HMCpar.NoMetropolisUntil =10;
HMCpar.MetropolisTest = true; HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed); sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed); pRNG.SeedFixedIntegers(ParSeed);
SU3::HotConfiguration(pRNG, U); SU3::HotConfiguration(pRNG, U);
} else if ( StartType == ColdStart ) { } else if ( StartType == ColdStart ) {
// Cold start // Cold start
HMCpar.NoMetropolisUntil =10;
HMCpar.MetropolisTest = true; HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed); sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed); pRNG.SeedFixedIntegers(ParSeed);
SU3::ColdConfiguration(pRNG, U); SU3::ColdConfiguration(pRNG, U);
} else if ( StartType == TepidStart ) { } else if ( StartType == TepidStart ) {
// Tepid start // Tepid start
HMCpar.NoMetropolisUntil =10;
HMCpar.MetropolisTest = true; HMCpar.MetropolisTest = true;
sRNG.SeedFixedIntegers(SerSeed); sRNG.SeedFixedIntegers(SerSeed);
pRNG.SeedFixedIntegers(ParSeed); pRNG.SeedFixedIntegers(ParSeed);
SU3::TepidConfiguration(pRNG, U); SU3::TepidConfiguration(pRNG, U);
} else if ( StartType == CheckpointStart ) { } else if ( StartType == CheckpointStart ) {
HMCpar.NoMetropolisUntil =10;
HMCpar.MetropolisTest = true; HMCpar.MetropolisTest = true;
// CheckpointRestart // CheckpointRestart
Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG); Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);

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

@@ -1,6 +1,6 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/utils/WilsonLoops.h Source file: ./lib/qcd/utils/WilsonLoops.h
@@ -25,472 +25,501 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
with this program; if not, write to the Free Software Foundation, Inc., with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#ifndef QCD_UTILS_WILSON_LOOPS_H #ifndef QCD_UTILS_WILSON_LOOPS_H
#define QCD_UTILS_WILSON_LOOPS_H #define QCD_UTILS_WILSON_LOOPS_H
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
// Common wilson loop observables
template<class Gimpl>
class WilsonLoops : public Gimpl {
public:
INHERIT_GIMPL_TYPES(Gimpl);
typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
//////////////////////////////////////////////////
// directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
// 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
//////////////////////////////////////////////////
static void traceDirPlaquette(LatticeComplex &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu)
{
GaugeMat sp(U[0]._grid);
dirPlaquette(sp,U,mu,nu);
plaq=trace(sp);
}
//////////////////////////////////////////////////
// sum over all planes of plaquette
//////////////////////////////////////////////////
static void sitePlaquette(LatticeComplex &Plaq,const std::vector<GaugeMat> &U)
{
LatticeComplex sitePlaq(U[0]._grid);
Plaq=zero;
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
traceDirPlaquette(sitePlaq,U,mu,nu);
Plaq = Plaq + sitePlaq;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumPlaquette(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
LatticeComplex Plaq(Umu._grid);
sitePlaquette(Plaq,U);
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgPlaquette(const GaugeLorentz &Umu){
RealD sumplaq = sumPlaquette(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0*Nd*(Nd-1))/2.0;
return sumplaq/vol/faces/Nc; // Nd , Nc dependent... FIXME
}
////////////////////////////////////////////////// // Common wilson loop observables
// average over traced single links template <class Gimpl> class WilsonLoops : public Gimpl {
////////////////////////////////////////////////// public:
static RealD linkTrace(const GaugeLorentz &Umu){ INHERIT_GIMPL_TYPES(Gimpl);
std::vector<GaugeMat> U(4,Umu._grid);
LatticeComplex Tr(Umu._grid); Tr=zero;
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
Tr = Tr+trace(U[mu]);
}
TComplex Tp = sum(Tr);
Complex p = TensorRemove(Tp);
double vol = Umu._grid->gSites();
return p.real()/vol/4.0/3.0;
};
//////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple,const GaugeLorentz &Umu,int mu, int nu){
GridBase *grid = Umu._grid; typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
std::vector<GaugeMat> U(4,grid);
for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d);
}
staple = zero;
if(nu != mu) {
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple+=Gimpl::ShiftStaple( //////////////////////////////////////////////////
Gimpl::CovShiftForward (U[nu],nu, // directed plaquette oriented in mu,nu plane
Gimpl::CovShiftBackward(U[mu],mu, //////////////////////////////////////////////////
Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu); static void dirPlaquette(GaugeMat &plaq, const std::vector<GaugeMat> &U,
const int mu, const int nu) {
// __ // 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
staple+=Gimpl::ShiftStaple( // purpose of deriving
Gimpl::CovShiftBackward(U[nu],nu, // from Gimpl.
Gimpl::CovShiftBackward(U[mu],mu,U[nu])),mu); plaq = Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftBackward(
} U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu])));
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple,const GaugeLorentz &Umu,int mu){
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd,grid);
for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d);
}
staple = zero;
GaugeMat tmp(grid);
for(int nu=0;nu<Nd;nu++){
if(nu != mu) {
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple+=Gimpl::ShiftStaple(
Gimpl::CovShiftForward (U[nu],nu,
Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu);
// __
// |
// |__
//
//
staple+=Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu],nu,
Gimpl::CovShiftBackward(U[mu],mu,U[nu])),mu);
}
}
} }
//////////////////////////////////////////////////
// trace of directed plaquette oriented in mu,nu plane
////////////////////////////////////////////////// //////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu, upper part static void traceDirPlaquette(LatticeComplex &plaq,
////////////////////////////////////////////////// const std::vector<GaugeMat> &U, const int mu,
static void StapleUpper(GaugeMat &staple,const GaugeLorentz &Umu,int mu, int nu){ const int nu) {
staple = zero;
if(nu != mu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4,grid);
for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d);
}
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple+=Gimpl::ShiftStaple(
Gimpl::CovShiftForward (U[nu],nu,
Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftIdentityBackward(U[nu],nu))),mu);
}
}
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////
static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
rect = Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(Gimpl::CovShiftForward(U[nu],nu,Gimpl::CovShiftForward(U[mu],mu,U[mu]))) ;
rect = rect +
Gimpl::CovShiftForward(U[mu],mu,Gimpl::CovShiftForward(U[nu],nu,U[nu]))* // ->||
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)
{
GaugeMat sp(U[0]._grid); GaugeMat sp(U[0]._grid);
dirRectangle(sp,U,mu,nu); dirPlaquette(sp, U, mu, nu);
rect=trace(sp); plaq = trace(sp);
} }
static void siteRectangle(LatticeComplex &Rect,const std::vector<GaugeMat> &U) //////////////////////////////////////////////////
{ // sum over all planes of plaquette
LatticeComplex siteRect(U[0]._grid); //////////////////////////////////////////////////
Rect=zero; static void sitePlaquette(LatticeComplex &Plaq,
for(int mu=1;mu<Nd;mu++){ const std::vector<GaugeMat> &U) {
for(int nu=0;nu<mu;nu++){ LatticeComplex sitePlaq(U[0]._grid);
traceDirRectangle(siteRect,U,mu,nu); Plaq = zero;
Rect = Rect + siteRect; for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) {
traceDirPlaquette(sitePlaq, U, mu, nu);
Plaq = Plaq + sitePlaq;
} }
} }
} }
//////////////////////////////////////////////////
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette // sum over all x,y,z,t and over all planes of plaquette
////////////////////////////////////////////////// //////////////////////////////////////////////////
static RealD sumRectangle(const GaugeLorentz &Umu){ static RealD sumPlaquette(const GaugeLorentz &Umu) {
std::vector<GaugeMat> U(Nd,Umu._grid); std::vector<GaugeMat> U(4, Umu._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);
} }
LatticeComplex Rect(Umu._grid); LatticeComplex Plaq(Umu._grid);
siteRectangle(Rect,U); sitePlaquette(Plaq, U);
TComplex Tp = sum(Rect); TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp); Complex p = TensorRemove(Tp);
return p.real(); return p.real();
} }
////////////////////////////////////////////////// //////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette // average over all x,y,z,t and over all planes of plaquette
////////////////////////////////////////////////// //////////////////////////////////////////////////
static RealD avgRectangle(const GaugeLorentz &Umu){ static RealD avgPlaquette(const GaugeLorentz &Umu) {
RealD sumplaq = sumPlaquette(Umu);
RealD sumrect = sumRectangle(Umu);
double vol = Umu._grid->gSites(); double vol = Umu._grid->gSites();
double faces = (1.0 * Nd * (Nd - 1)) / 2.0;
double faces = (1.0*Nd*(Nd-1)); // 2 distinct orientations summed return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME
}
return sumrect/vol/faces/Nc; // Nd , Nc dependent... FIXME
//////////////////////////////////////////////////
// average over traced single links
//////////////////////////////////////////////////
static RealD linkTrace(const GaugeLorentz &Umu) {
std::vector<GaugeMat> U(4, Umu._grid);
LatticeComplex Tr(Umu._grid);
Tr = zero;
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
Tr = Tr + trace(U[mu]);
}
TComplex Tp = sum(Tr);
Complex p = TensorRemove(Tp);
double vol = Umu._grid->gSites();
return p.real() / vol / 4.0 / 3.0;
};
//////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu,
int nu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
staple = zero;
if (nu != mu) {
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftForward(
U[nu], nu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
mu);
// __
// |
// |__
//
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
mu);
}
} }
////////////////////////////////////////////////// //////////////////////////////////////////////////
// the sum over all staples on each site // the sum over all staples on each site
////////////////////////////////////////////////// //////////////////////////////////////////////////
static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){ static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
U2 = U * Cshift(U,mu,1);
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
staple = zero;
GaugeMat tmp(grid);
for (int nu = 0; nu < Nd; nu++) {
if (nu != mu) {
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftForward(
U[nu], nu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
mu);
// __
// |
// |__
//
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
mu);
}
}
}
//////////////////////////////////////////////////
// the sum over all staples on each site in direction mu,nu, upper part
//////////////////////////////////////////////////
static void StapleUpper(GaugeMat &staple, const GaugeLorentz &Umu, int mu,
int nu) {
staple = zero;
if (nu != mu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
// mu
// ^
// |__> nu
// __
// |
// __|
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftForward(
U[nu], nu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
mu);
}
}
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////
static void dirRectangle(GaugeMat &rect, const std::vector<GaugeMat> &U,
const int mu, const int nu) {
rect = Gimpl::CovShiftForward(
U[mu], mu, Gimpl::CovShiftForward(U[mu], mu, U[nu])) * // ->->|
adj(Gimpl::CovShiftForward(
U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[mu])));
rect = rect +
Gimpl::CovShiftForward(
U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[nu])) * // ->||
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) {
GaugeMat sp(U[0]._grid);
dirRectangle(sp, U, mu, nu);
rect = trace(sp);
}
static void siteRectangle(LatticeComplex &Rect,
const std::vector<GaugeMat> &U) {
LatticeComplex siteRect(U[0]._grid);
Rect = zero;
for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) {
traceDirRectangle(siteRect, U, mu, nu);
Rect = Rect + siteRect;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumRectangle(const GaugeLorentz &Umu) {
std::vector<GaugeMat> U(Nd, Umu._grid);
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
}
LatticeComplex Rect(Umu._grid);
siteRectangle(Rect, U);
TComplex Tp = sum(Rect);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgRectangle(const GaugeLorentz &Umu) {
RealD sumrect = sumRectangle(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0 * Nd * (Nd - 1)); // 2 distinct orientations summed
return sumrect / vol / faces / Nc; // Nd , Nc dependent... FIXME
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void RectStapleDouble(GaugeMat &U2, const GaugeMat &U, int mu) {
U2 = U * Cshift(U, mu, 1);
} }
//////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////
// Hop by two optimisation strategy does not work nicely with Gparity. (could do, // 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). // 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 . // 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;
GridBase *grid = U[0]._grid; GridBase *grid = U[0]._grid;
GaugeMat Staple2x1 (grid); GaugeMat Staple2x1(grid);
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) {
// Up staple ___ ___ // Up staple ___ ___
// | | // | |
tmp = Cshift(adj(U[nu]),nu,-1); tmp = Cshift(adj(U[nu]), nu, -1);
tmp = adj(U2[mu])*tmp; tmp = adj(U2[mu]) * tmp;
tmp = Cshift(tmp,mu,-2); tmp = Cshift(tmp, mu, -2);
Staple2x1 = Gimpl::CovShiftForward (U[nu],nu,tmp); Staple2x1 = Gimpl::CovShiftForward(U[nu], nu, tmp);
// Down staple
// |___ ___|
//
tmp = adj(U2[mu]) * U[nu];
Staple2x1 += Gimpl::CovShiftBackward(U[nu], nu, Cshift(tmp, mu, -2));
// Down staple // ___ ___
// |___ ___| // | ___|
// // |___ ___|
tmp = adj(U2[mu])*U[nu]; //
Staple2x1+= Gimpl::CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
Stap += Cshift(Gimpl::CovShiftForward(U[mu], mu, Staple2x1), mu, 1);
// ___ ___ // ___ ___
// | ___| // |___ |
// |___ ___| // |___ ___|
// //
Stap+= Cshift(Gimpl::CovShiftForward (U[mu],mu,Staple2x1),mu,1); // tmp= Staple2x1* Cshift(U[mu],mu,-2);
// Stap+= Cshift(tmp,mu,1) ;
Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1);
;
// ___ ___ // --
// |___ | // | |
// |___ ___| //
// // | |
// tmp= Staple2x1* Cshift(U[mu],mu,-2); tmp = Cshift(adj(U2[nu]), nu, -2);
// Stap+= Cshift(tmp,mu,1) ; tmp = Gimpl::CovShiftBackward(U[mu], mu, tmp);
Stap+= Cshift(Staple2x1,mu,1)*Cshift(U[mu],mu,-1); ; tmp = U2[nu] * Cshift(tmp, nu, 2);
Stap += Cshift(tmp, mu, 1);
// -- // | |
// | | //
// // | |
// | | // --
tmp = Cshift(adj(U2[nu]),nu,-2);
tmp = Gimpl::CovShiftBackward(U[mu],mu,tmp);
tmp = U2[nu]*Cshift(tmp,nu,2);
Stap+= Cshift(tmp, mu, 1);
// | | tmp = Gimpl::CovShiftBackward(U[mu], mu, U2[nu]);
// tmp = adj(U2[nu]) * tmp;
// | | tmp = Cshift(tmp, nu, -2);
// -- Stap += Cshift(tmp, mu, 1);
}
tmp = Gimpl::CovShiftBackward(U[mu],mu,U2[nu]);
tmp = adj(U2[nu])*tmp;
tmp = Cshift(tmp,nu,-2);
Stap+=Cshift(tmp, mu, 1);
}}
}
static void RectStaple(GaugeMat &Stap,const GaugeLorentz & Umu,int mu)
{
RectStapleUnoptimised(Stap,Umu,mu);
}
static void RectStaple(const GaugeLorentz & Umu,GaugeMat &Stap,
std::vector<GaugeMat> &U2,
std::vector<GaugeMat> &U, int mu)
{
if ( Gimpl::isPeriodicGaugeField() ){
RectStapleOptimised(Stap,U2,U,mu);
} else {
RectStapleUnoptimised(Stap,Umu,mu);
} }
} }
static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){ static void RectStaple(GaugeMat &Stap, const GaugeLorentz &Umu, int mu) {
RectStapleUnoptimised(Stap, Umu, mu);
}
static void RectStaple(const GaugeLorentz &Umu, GaugeMat &Stap,
std::vector<GaugeMat> &U2, std::vector<GaugeMat> &U,
int mu) {
if (Gimpl::isPeriodicGaugeField()) {
RectStapleOptimised(Stap, U2, U, mu);
} else {
RectStapleUnoptimised(Stap, Umu, mu);
}
}
static void RectStapleUnoptimised(GaugeMat &Stap, const GaugeLorentz &Umu,
int mu) {
GridBase *grid = Umu._grid; GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd,grid); std::vector<GaugeMat> U(Nd, grid);
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);
} }
Stap=zero; Stap = zero;
for(int nu=0;nu<Nd;nu++){ for (int nu = 0; nu < Nd; nu++) {
if ( nu!=mu) { if (nu != mu) {
// __ ___ // __ ___
// | __ | // | __ |
// //
Stap+= Gimpl::ShiftStaple( Stap += Gimpl::ShiftStaple(
Gimpl::CovShiftForward (U[mu],mu, Gimpl::CovShiftForward(
Gimpl::CovShiftForward (U[nu],nu, U[mu], mu,
Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftForward(
Gimpl::CovShiftBackward(U[mu],mu, U[nu], nu,
Gimpl::CovShiftIdentityBackward(U[nu],nu))))) , mu); Gimpl::CovShiftBackward(
U[mu], mu,
Gimpl::CovShiftBackward(
U[mu], mu,
Gimpl::CovShiftIdentityBackward(U[nu], nu))))),
mu);
// __ // __
// |__ __ | // |__ __ |
Stap+= Gimpl::ShiftStaple( Stap += Gimpl::ShiftStaple(
Gimpl::CovShiftForward (U[mu],mu, Gimpl::CovShiftForward(
Gimpl::CovShiftBackward(U[nu],nu, U[mu], mu,
Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(
Gimpl::CovShiftBackward(U[mu],mu, U[nu])))) , mu); U[nu], nu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])))),
mu);
// __ // __
// |__ __ | // |__ __ |
Stap+= Gimpl::ShiftStaple( Stap += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu],nu, Gimpl::CovShiftBackward(
Gimpl::CovShiftBackward(U[mu],mu, U[nu], nu,
Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(
Gimpl::CovShiftForward(U[nu],nu,U[mu])))) , mu); U[mu], mu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[mu])))),
mu);
// __ ___ // __ ___
// |__ | // |__ |
Stap+= Gimpl::ShiftStaple( Stap += Gimpl::ShiftStaple(
Gimpl::CovShiftForward (U[nu],nu, Gimpl::CovShiftForward(
Gimpl::CovShiftBackward(U[mu],mu, U[nu], nu,
Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(
Gimpl::CovShiftBackward(U[nu],nu,U[mu])))) , mu); U[mu], mu,
Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftBackward(U[nu], nu, U[mu])))),
mu);
// -- // --
// | | // | |
// //
// | | // | |
Stap+= Gimpl::ShiftStaple(
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))))) , mu);
Stap += Gimpl::ShiftStaple(
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))))),
mu);
// | | // | |
// //
// | | // | |
// -- // --
Stap+= Gimpl::ShiftStaple( Stap += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu],nu, Gimpl::CovShiftBackward(
Gimpl::CovShiftBackward(U[nu],nu, U[nu], nu,
Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftBackward(
Gimpl::CovShiftForward (U[nu],nu,U[nu])))) , mu); U[nu], nu,
}} Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[nu])))),
mu);
}
}
} }
}; };
typedef WilsonLoops<PeriodicGimplR> ColourWilsonLoops;
typedef WilsonLoops<PeriodicGimplR> ColourWilsonLoops; typedef WilsonLoops<PeriodicGimplR> U1WilsonLoops;
typedef WilsonLoops<PeriodicGimplR> U1WilsonLoops; typedef WilsonLoops<PeriodicGimplR> SU2WilsonLoops;
typedef WilsonLoops<PeriodicGimplR> SU2WilsonLoops; typedef WilsonLoops<PeriodicGimplR> SU3WilsonLoops;
typedef WilsonLoops<PeriodicGimplR> SU3WilsonLoops; }
}
}}
#endif #endif