1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-14 01:35:36 +00:00

Merge branch 'feature/hadrons' into feature/rare_kaon

This commit is contained in:
Lanny91 2017-02-01 10:13:29 +00:00
commit 97053adcb5
84 changed files with 5028 additions and 1125 deletions

12
.gitignore vendored
View File

@ -105,3 +105,15 @@ lib/fftw/*
################## ##################
m4/lt* m4/lt*
m4/libtool.m4 m4/libtool.m4
# Buck files #
##############
.buck*
buck-out
BUCK
make-bin-BUCK.sh
# generated sources #
#####################
lib/qcd/spin/gamma-gen/*.h
lib/qcd/spin/gamma-gen/*.cc

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR; typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;
@ -321,7 +321,7 @@ int main (int argc, char ** argv)
ref = zero; ref = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
// ref = src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x // ref = src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
tmp = U[mu]*Cshift(src,mu+1,1); tmp = U[mu]*Cshift(src,mu+1,1);
for(int i=0;i<ref._odata.size();i++){ for(int i=0;i<ref._odata.size();i++){
ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ; ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
void benchDw(std::vector<int> & L, int Ls, int threads, int report =0 ); void benchDw(std::vector<int> & L, int Ls, int threads, int report =0 );

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
bool overlapComms = false; bool overlapComms = false;
@ -106,7 +106,7 @@ int main (int argc, char ** argv)
{ // Naive wilson implementation { // Naive wilson implementation
ref = zero; ref = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
// ref = src + Gamma(Gamma::GammaX)* src ; // 1-gamma_x // ref = src + Gamma(Gamma::Algebra::GammaX)* src ; // 1-gamma_x
tmp = U[mu]*Cshift(src,mu,1); tmp = U[mu]*Cshift(src,mu,1);
for(int i=0;i<ref._odata.size();i++){ for(int i=0;i<ref._odata.size();i++){
ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ; ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
@ -159,7 +159,7 @@ int main (int argc, char ** argv)
ref = zero; ref = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
// ref = src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x // ref = src - Gamma(Gamma::Algebra::GammaX)* src ; // 1+gamma_x
tmp = U[mu]*Cshift(src,mu,1); tmp = U[mu]*Cshift(src,mu,1);
for(int i=0;i<ref._odata.size();i++){ for(int i=0;i<ref._odata.size();i++){
ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ; ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;

View File

@ -30,11 +30,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
bool overlapComms = false; bool overlapComms = false;

View File

@ -319,7 +319,7 @@ AM_CONDITIONAL(BUILD_COMMS_MPI3L, [ test "${comms_type}X" == "mpi3lX" ] )
AM_CONDITIONAL(BUILD_COMMS_NONE, [ test "${comms_type}X" == "noneX" ]) AM_CONDITIONAL(BUILD_COMMS_NONE, [ test "${comms_type}X" == "noneX" ])
############### RNG selection ############### RNG selection
AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937],\ AC_ARG_ENABLE([rng],[AC_HELP_STRING([--enable-rng=ranlux48|mt19937|sitmo],\
[Select Random Number Generator to be used])],\ [Select Random Number Generator to be used])],\
[ac_RNG=${enable_rng}],[ac_RNG=ranlux48]) [ac_RNG=${enable_rng}],[ac_RNG=ranlux48])
@ -330,6 +330,9 @@ case ${ac_RNG} in
mt19937) mt19937)
AC_DEFINE([RNG_MT19937],[1],[RNG_MT19937] ) AC_DEFINE([RNG_MT19937],[1],[RNG_MT19937] )
;; ;;
sitmo)
AC_DEFINE([RNG_SITMO],[1],[RNG_SITMO] )
;;
*) *)
AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]); AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]);
;; ;;

View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MContraction/Meson.hpp
Copyright (C) 2015 Copyright (C) 2015
Copyright (C) 2016 Copyright (C) 2016
Copyright (C) 2017
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Andrew Lawson <andrew.lawson1991@gmail.com> Andrew Lawson <andrew.lawson1991@gmail.com>
@ -36,17 +37,17 @@ See the full license in the file "LICENSE" in the top level distribution directo
#include <Grid/Hadrons/ModuleFactory.hpp> #include <Grid/Hadrons/ModuleFactory.hpp>
namespace Grid { namespace Grid {
// Overload >> to extract gamma pair from "[g1 g2]" string. // Overload >> to extract gamma pair from "<g1 g2>" string.
template <typename T1, typename T2> template <typename T1, typename T2>
inline std::istringstream &operator>>(std::istringstream &sstr, inline std::istringstream &operator>>(std::istringstream &sstr,
std::pair<T1, T2> &buf) std::pair<T1, T2> &buf)
{ {
T1 buf1; unsigned int buf1;
T2 buf2; unsigned int buf2;
char c; char c;
sstr >> c >> buf1 >> buf2 >> c; sstr >> c >> buf1 >> buf2 >> c;
sstr.peek(); sstr.peek();
buf = std::make_pair(buf1, buf2); buf = std::make_pair((T1)buf1, (T2)buf2);
return sstr; return sstr;
} }
} }
@ -66,15 +67,16 @@ BEGIN_HADRONS_NAMESPACE
in a sequence (e.g. "[15 7][7 15][7 7]"). in a sequence (e.g. "[15 7][7 15][7 7]").
Special values: "all" - perform all possible contractions. Special values: "all" - perform all possible contractions.
- mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."),
*/ given as multiples of (2*pi) / L.
*/
/****************************************************************************** /******************************************************************************
* TMeson * * TMeson *
******************************************************************************/ ******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction) BEGIN_MODULE_NAMESPACE(MContraction)
typedef std::pair<unsigned int, unsigned int> GammaPair; typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
class MesonPar: Serializable class MesonPar: Serializable
{ {
@ -83,6 +85,7 @@ public:
std::string, q1, std::string, q1,
std::string, q2, std::string, q2,
std::string, gammas, std::string, gammas,
std::string, mom,
std::string, output); std::string, output);
}; };
@ -96,8 +99,8 @@ public:
{ {
public: public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result, GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
unsigned int, gamma_snk, Gamma::Algebra, gamma_snk,
unsigned int, gamma_src, Gamma::Algebra, gamma_src,
std::vector<Complex>, corr); std::vector<Complex>, corr);
}; };
public: public:
@ -148,13 +151,14 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
if (par().gammas.compare("all") == 0) if (par().gammas.compare("all") == 0)
{ {
// Do all contractions. // Do all contractions.
unsigned int n_gam = Ns*Ns; unsigned int n_gam = Ns * Ns;
gammaList.resize(n_gam*n_gam); gammaList.resize(n_gam*n_gam);
for (unsigned int i = 0; i < n_gam; ++i) for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
{ {
for (unsigned int j = 0; j < n_gam; ++j) for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
{ {
gammaList.push_back(std::make_pair(i, j)); gammaList.push_back(std::make_pair((Gamma::Algebra)i,
(Gamma::Algebra)j));
} }
} }
} }
@ -178,22 +182,31 @@ void TMeson<FImpl1, FImpl2>::execute(void)
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1); PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2); PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
LatticeComplex c(env().getGrid()); LatticeComplex c(env().getGrid());
SpinMatrix g[Ns*Ns], g5; Gamma g5(Gamma::Algebra::Gamma5);
std::vector<GammaPair> gammaList; std::vector<GammaPair> gammaList;
std::vector<TComplex> buf; std::vector<TComplex> buf;
std::vector<Result> result; std::vector<Result> result;
std::vector<Real> p;
g5 = makeGammaProd(Ns*Ns - 1); p = strToVec<Real>(par().mom);
for (int i = 0; i < Ns*Ns; ++i) LatticeComplex ph(env().getGrid()), coor(env().getGrid());
Complex i(0.0,1.0);
ph = zero;
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{ {
g[i] = makeGammaProd(i); LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
} }
ph = exp(-2*M_PI*i*ph);
parseGammaString(gammaList); parseGammaString(gammaList);
result.resize(gammaList.size()); result.resize(gammaList.size());
for (unsigned int i = 0; i < result.size(); ++i) for (unsigned int i = 0; i < result.size(); ++i)
{ {
c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5); Gamma gSnk(gammaList[i].first);
Gamma gSrc(gammaList[i].second);
c = trace((g5*gSnk)*q1*(gSrc*g5)*adj(q2))*ph;
sliceSum(c, buf, Tp); sliceSum(c, buf, Tp);
result[i].gamma_snk = gammaList[i].first; result[i].gamma_snk = gammaList[i].first;

View File

@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MSink/Wall.hpp Source file: extras/Hadrons/Modules/MSink/Wall.hpp
Copyright (C) 2016 Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com> Author: Andrew Lawson <andrew.lawson1991@gmail.com>
@ -131,9 +131,9 @@ void TWall<FImpl>::execute(void)
for(unsigned int mu = 0; mu < Nd; mu++) for(unsigned int mu = 0; mu < Nd; mu++)
{ {
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor; ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
} }
ph = exp(-i*ph); ph = exp(-2*M_PI*i*ph);
sliceSum<SitePropagator>(ph*q, prop, Tp); sliceSum<SitePropagator>(ph*q, prop, Tp);
} }

View File

@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp
Copyright (C) 2015 Copyright (C) 2015
Copyright (C) 2016 Copyright (C) 2016
Copyright (C) 2017
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
@ -63,7 +64,7 @@ public:
std::string, q, std::string, q,
unsigned int, tA, unsigned int, tA,
unsigned int, tB, unsigned int, tB,
unsigned int, gamma, Gamma::Algebra, gamma,
std::string, mom); std::string, mom);
}; };
@ -140,21 +141,20 @@ void TSeqGamma<FImpl>::execute(void)
PropagatorField &q = *env().template getObject<PropagatorField>(par().q); PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
Lattice<iScalar<vInteger>> t(env().getGrid()); Lattice<iScalar<vInteger>> t(env().getGrid());
LatticeComplex ph(env().getGrid()), coor(env().getGrid()); LatticeComplex ph(env().getGrid()), coor(env().getGrid());
SpinMatrix g; Gamma g(par().gamma);
std::vector<Real> p; std::vector<Real> p;
Complex i(0.0,1.0); Complex i(0.0,1.0);
g = makeGammaProd(par().gamma);
p = strToVec<Real>(par().mom); p = strToVec<Real>(par().mom);
ph = zero; ph = zero;
for(unsigned int mu = 0; mu < env().getNd(); mu++) for(unsigned int mu = 0; mu < env().getNd(); mu++)
{ {
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor; ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
} }
ph = exp(i*ph); ph = exp(2*M_PI*i*ph);
LatticeCoordinate(t, Tp); LatticeCoordinate(t, Tp);
src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q); src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
} }
END_MODULE_NAMESPACE END_MODULE_NAMESPACE

View File

@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
Source file: extras/Hadrons/Modules/MSource/Wall.hpp Source file: extras/Hadrons/Modules/MSource/Wall.hpp
Copyright (C) 2016 Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com> Author: Andrew Lawson <andrew.lawson1991@gmail.com>
@ -132,9 +132,9 @@ void TWall<FImpl>::execute(void)
for(unsigned int mu = 0; mu < Nd; mu++) for(unsigned int mu = 0; mu < Nd; mu++)
{ {
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
ph = ph + p[mu]*coor; ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
} }
ph = exp(i*ph); ph = exp(2*M_PI*i*ph);
LatticeCoordinate(t, Tp); LatticeCoordinate(t, Tp);
src = 1.; src = 1.;
src = where((t == par().tW), src*ph, 0.*src); src = where((t == par().tW), src*ph, 0.*src);

View File

@ -41,7 +41,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#include <signal.h> #include <signal.h>
#include <iostream> #include <iostream>
#include <iterator> #include <iterator>
#include <Grid.h> #include <Grid/Grid.h>
#include <algorithm> #include <algorithm>
#include <iterator> #include <iterator>
#include <cstdlib> #include <cstdlib>

View File

@ -29,7 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
#include <cxxabi.h> #include <cxxabi.h>

View File

@ -26,8 +26,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
#include <PerfCount.h> #include <Grid/PerfCount.h>
namespace Grid { namespace Grid {

View File

@ -1,6 +1,6 @@
#include <Grid.h> #include <Grid/Grid.h>
#include <PerfCount.h> #include <Grid/PerfCount.h>
#include <Stat.h> #include <Grid/Stat.h>
namespace Grid { namespace Grid {

View File

@ -25,7 +25,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
double MultiShiftFunction::approx(double x) double MultiShiftFunction::approx(double x)

View File

@ -20,7 +20,7 @@
#include<iomanip> #include<iomanip>
#include<cassert> #include<cassert>
#include<algorithms/approx/Remez.h> #include<Grid/algorithms/approx/Remez.h>
// Constructor // Constructor
AlgRemez::AlgRemez(double lower, double upper, long precision) AlgRemez::AlgRemez(double lower, double upper, long precision)

View File

@ -36,7 +36,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
#include <iomanip> #include <iomanip>
#include <complex> #include <complex>
#include <typeinfo> #include <typeinfo>
#include <Grid.h> #include <Grid/Grid.h>
/** Sign function **/ /** Sign function **/

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#include <mpi.h> #include <mpi.h>
namespace Grid { namespace Grid {

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#include <mpi.h> #include <mpi.h>
namespace Grid { namespace Grid {

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
/////////////////////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////////////////////

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#include <mpp/shmem.h> #include <mpp/shmem.h>
namespace Grid { namespace Grid {

View File

@ -30,6 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
#define GRID_LATTICE_RNG_H #define GRID_LATTICE_RNG_H
#include <random> #include <random>
#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
namespace Grid { namespace Grid {
@ -114,10 +115,14 @@ namespace Grid {
typedef uint64_t RngStateType; typedef uint64_t RngStateType;
typedef std::ranlux48 RngEngine; typedef std::ranlux48 RngEngine;
static const int RngStateCount = 15; static const int RngStateCount = 15;
#else #elif RNG_MT19937
typedef std::mt19937 RngEngine; typedef std::mt19937 RngEngine;
typedef uint32_t RngStateType; typedef uint32_t RngStateType;
static const int RngStateCount = std::mt19937::state_size; static const int RngStateCount = std::mt19937::state_size;
#elif RNG_SITMO
typedef sitmo::prng_engine RngEngine;
typedef uint64_t RngStateType;
static const int RngStateCount = 4;
#endif #endif
std::vector<RngEngine> _generators; std::vector<RngEngine> _generators;
std::vector<std::uniform_real_distribution<RealD>> _uniform; std::vector<std::uniform_real_distribution<RealD>> _uniform;

View File

@ -14,7 +14,7 @@
#ifndef SOURCE_PUGIXML_CPP #ifndef SOURCE_PUGIXML_CPP
#define SOURCE_PUGIXML_CPP #define SOURCE_PUGIXML_CPP
#include <pugixml/pugixml.h> #include <Grid/pugixml/pugixml.h>
#include <stdlib.h> #include <stdlib.h>
#include <stdio.h> #include <stdio.h>

View File

@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#include <Grid/Eigen/Dense> #include <Grid/Eigen/Dense>
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -29,7 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#include <Grid/Eigen/Dense> #include <Grid/Eigen/Dense>
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -29,7 +29,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -30,7 +30,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -26,7 +26,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {

View File

@ -26,7 +26,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {

View File

@ -29,7 +29,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
@ -149,11 +149,11 @@ void WilsonFermion<Impl>::MeooeDag(const FermionField &in, FermionField &out) {
typedef Lattice<iSinglet<vector_type> > LatComplex; typedef Lattice<iSinglet<vector_type> > LatComplex;
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
std::vector<int> latt_size = _grid->_fdimensions; std::vector<int> latt_size = _grid->_fdimensions;

View File

@ -29,8 +29,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#include <PerfCount.h> #include <Grid/PerfCount.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
@ -503,11 +503,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
typedef iSinglet<ScalComplex> Tcomplex; typedef iSinglet<ScalComplex> Tcomplex;
typedef Lattice<iSinglet<vector_type> > LatComplex; typedef Lattice<iSinglet<vector_type> > LatComplex;
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
std::vector<int> latt_size = _grid->_fdimensions; std::vector<int> latt_size = _grid->_fdimensions;
@ -574,11 +574,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass)
{ {
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
GridBase *_grid = _FourDimGrid; GridBase *_grid = _FourDimGrid;

View File

@ -28,7 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
directory directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {

View File

@ -30,7 +30,7 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid.h> #include <Grid/Grid.h>
namespace Grid { namespace Grid {

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#define REGISTER #define REGISTER

View File

@ -25,7 +25,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {

View File

@ -80,7 +80,7 @@ class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat; Matrix &_Mat;
Gamma g5; Gamma g5;
public: public:
Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Gamma5) {}; Gamma5HermitianLinearOperator(Matrix &Mat): _Mat(Mat), g5(Gamma::Algebra::Gamma5) {};
void Op (const Field &in, Field &out){ void Op (const Field &in, Field &out){
HermOp(in,out); HermOp(in,out);
} }

View File

@ -27,7 +27,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid{ namespace Grid{
namespace QCD{ namespace QCD{

View File

@ -1,95 +1,72 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/spin/Dirac.cc Source file: lib/qcd/spin/Dirac.cc
Copyright (C) 2015 Copyright (C) 2015
Copyright (C) 2016
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <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/Grid.h>
namespace Grid { namespace Grid {
namespace QCD {
namespace QCD { #include "GammaMulTable.h"
Gamma::GammaMatrix Gamma::GammaMatrices [] = { const std::array<const char *, Gamma::nGamma> Gamma::name = {{
Gamma::Identity,
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT,
Gamma::Gamma5,
Gamma::MinusIdentity,
Gamma::MinusGammaX,
Gamma::MinusGammaY,
Gamma::MinusGammaZ,
Gamma::MinusGammaT,
Gamma::MinusGamma5
};
const char *Gamma::GammaMatrixNames[] = {
"Identity ",
"GammaX ",
"GammaY ",
"GammaZ ",
"GammaT ",
"Gamma5 ",
"-Identity",
"-GammaX ",
"-GammaY ",
"-GammaZ ",
"-GammaT ",
"-Gamma5 ", "-Gamma5 ",
" " "Gamma5 ",
}; "-GammaT ",
"GammaT ",
"-GammaTGamma5",
"GammaTGamma5 ",
"-GammaX ",
"GammaX ",
"-GammaXGamma5",
"GammaXGamma5 ",
"-GammaY ",
"GammaY ",
"-GammaYGamma5",
"GammaYGamma5 ",
"-GammaZ ",
"GammaZ ",
"-GammaZGamma5",
"GammaZGamma5 ",
"-Identity ",
"Identity ",
"-SigmaXT ",
"SigmaXT ",
"-SigmaXY ",
"SigmaXY ",
"-SigmaXZ ",
"SigmaXZ ",
"-SigmaYT ",
"SigmaYT ",
"-SigmaYZ ",
"SigmaYZ ",
"-SigmaZT ",
"SigmaZT "}};
SpinMatrix makeGammaProd(const unsigned int i) }}
{
SpinMatrix g;
g = 1.;
if (i & 0x1)
{
g = g*Gamma(Gamma::GammaMatrix::GammaX);
}
if (i & 0x2)
{
g = g*Gamma(Gamma::GammaMatrix::GammaY);
}
if (i & 0x4)
{
g = g*Gamma(Gamma::GammaMatrix::GammaZ);
}
if (i & 0x8)
{
g = g*Gamma(Gamma::GammaMatrix::GammaT);
}
return g;
}
// void sprojMul( vHalfSpinColourVector &out,vColourMatrix &u, vSpinColourVector &in){
// vHalfSpinColourVector hspin;
// spProjXp(hspin,in);
// mult(&out,&u,&hspin);
// }
}
}

File diff suppressed because it is too large Load Diff

1121
lib/qcd/spin/GammaMulTable.h Normal file

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -47,7 +47,7 @@ void axpibg5x(Lattice<vobj> &z,const Lattice<vobj> &x,Coeff a,Coeff b)
GridBase *grid=x._grid; GridBase *grid=x._grid;
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss++){ for(int ss=0;ss<grid->oSites();ss++){
vobj tmp; vobj tmp;
@ -80,7 +80,7 @@ void ag5xpby_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
conformable(x,z); conformable(x,z);
GridBase *grid=x._grid; GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0]; int Ls = grid->_rdimensions[0];
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp; vobj tmp;
@ -98,7 +98,7 @@ void axpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const L
conformable(x,z); conformable(x,z);
GridBase *grid=x._grid; GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0]; int Ls = grid->_rdimensions[0];
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp; vobj tmp;
@ -116,7 +116,7 @@ void ag5xpbg5y_ssp(Lattice<vobj> &z,Coeff a,const Lattice<vobj> &x,Coeff b,const
conformable(x,z); conformable(x,z);
GridBase *grid=x._grid; GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0]; int Ls = grid->_rdimensions[0];
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp1; vobj tmp1;
@ -168,7 +168,7 @@ void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
z.checkerboard = x.checkerboard; z.checkerboard = x.checkerboard;
conformable(x,z); conformable(x,z);
int Ls = grid->_rdimensions[0]; int Ls = grid->_rdimensions[0];
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp; vobj tmp;

View File

@ -25,7 +25,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
using namespace Grid; using namespace Grid;
using namespace std; using namespace std;

View File

@ -1,4 +1,4 @@
#include <Grid.h> #include <Grid/Grid.h>
using namespace Grid; using namespace Grid;
#ifndef H5_NO_NAMESPACE #ifndef H5_NO_NAMESPACE

View File

@ -150,14 +150,14 @@ friend inline bool operator==(const cname &lhs, const cname &rhs) {\
class name: public Grid::Serializable\ class name: public Grid::Serializable\
{\ {\
public:\ public:\
enum EnumType\ enum\
{\ {\
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\ GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
undefname = -1\ undefname = -1\
};\ };\
public:\ public:\
name(void): value_(undefname) {};\ name(void): value_(undefname) {};\
name(EnumType value): value_(value) {};\ name(int value): value_(value) {};\
template <typename T>\ template <typename T>\
static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\ static inline void write(Grid::Writer<T> &WR,const std::string &s, const name &obj)\
{\ {\
@ -177,7 +177,7 @@ public:\
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\ GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
else {obj = name::undefname;}\ else {obj = name::undefname;}\
}\ }\
inline operator EnumType(void) const\ inline operator int(void) const\
{\ {\
return value_;\ return value_;\
}\ }\
@ -190,7 +190,7 @@ public:\
return os;\ return os;\
}\ }\
private:\ private:\
EnumType value_;\ int value_;\
}; };

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
using namespace Grid; using namespace Grid;
using namespace std; using namespace std;

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
using namespace Grid; using namespace Grid;
using namespace std; using namespace std;

View File

@ -0,0 +1,390 @@
// Copyright (c) 2012-2016 M.A. (Thijs) van den Berg, http://sitmo.com/
//
// Use, modification and distribution are subject to the MIT Software License.
//
// The MIT License (MIT)
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.
// version history:
// version 1, 6 Sep 2012
// version 2, 10 Dec 2013
// bug fix in the discard() routine, it was discarding to many elements
// added the version() method
// version 3...5, 13 Dec 2013
// fixed type-conversion earning
// fixed potential issues with constructor template matching
// version 6, 4 March 2016
// made min() max() constexpr for C+11 compiler (thanks to James Joseph Balamuta)
#ifndef SITMO_PRNG_ENGINE_HPP
#define SITMO_PRNG_ENGINE_HPP
#include <iostream>
#ifdef __GNUC__
#include <stdint.h> // respecting the C99 standard.
#endif
#ifdef _MSC_VER
typedef unsigned __int64 uint64_t; // Visual Studio 6.0(VC6) and newer..
typedef unsigned __int32 uint32_t;
#endif
// Double mixing function
#define MIX2(x0,x1,rx,z0,z1,rz) \
x0 += x1; \
z0 += z1; \
x1 = (x1 << rx) | (x1 >> (64-rx)); \
z1 = (z1 << rz) | (z1 >> (64-rz)); \
x1 ^= x0; \
z1 ^= z0;
// Double mixing function with key adition
#define MIXK(x0,x1,rx,z0,z1,rz,k0,k1,l0,l1) \
x1 += k1; \
z1 += l1; \
x0 += x1+k0; \
z0 += z1+l0; \
x1 = (x1 << rx) | (x1 >> (64-rx)); \
z1 = (z1 << rz) | (z1 >> (64-rz)); \
x1 ^= x0; \
z1 ^= z0; \
namespace sitmo {
// enable_if for C__98 compilers
template<bool C, typename T = void>
struct sitmo_enable_if { typedef T type; };
template<typename T>
struct sitmo_enable_if<false, T> { };
// SFINAE check for the existence of a "void generate(int*,int*)"member function
template<typename T>
struct has_generate_template
{
typedef char (&Two)[2];;
template<typename F, void (F::*)(int *, int *)> struct helper {};
template<typename C> static char test(helper<C, &C::template generate<int*> >*);
template<typename C> static Two test(...);
static bool const value = sizeof(test<T>(0)) == sizeof(char);
};
class prng_engine
{
public:
// "req" are requirements as stated in the C++ 11 draft n3242=11-0012
//
// req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 1
typedef uint32_t result_type;
// req: 26.5.1.3 Uniform random number generator requirements, p.906, table 116, row 3 & 4
#if __cplusplus <= 199711L
static result_type (min)() { return 0; }
static result_type (max)() { return 0xFFFFFFFF; }
#else
static constexpr result_type (min)() { return 0; }
static constexpr result_type (max)() { return 0xFFFFFFFF; }
#endif
// -------------------------------------------------
// Constructors
// -------------------------------------------------
// req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 1
// Creates an engine with the same initial state as all other
// default-constructed engines of type E.
prng_engine()
{
seed();
}
// req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 2
// Creates an engine that compares equal to x.
prng_engine(const prng_engine& x)
{
for (unsigned short i=0; i<4; ++i) {
_s[i] = x._s[i];
_k[i] = x._k[i];
_o[i] = x._o[i];
}
_o_counter = x._o_counter;
}
// req: 26.5.1.4 Random number engine requirements, p.907 table 117, row 3
// Creates an engine with initial O(size of state) state determined by s.
prng_engine(uint32_t s)
{
seed(s);
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 4
// Creates an engine with an initial state that depends on a sequence
// produced by one call to q.generate.
template<class Seq>
prng_engine(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
{
seed(q);
}
// -------------------------------------------------
// Seeding
// -------------------------------------------------
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 5
void seed()
{
for (unsigned short i=0; i<4; ++i) {
_k[i] = 0;
_s[i] = 0;
}
_o_counter = 0;
_o[0] = 0x09218ebde6c85537;
_o[1] = 0x55941f5266d86105;
_o[2] = 0x4bd25e16282434dc;
_o[3] = 0xee29ec846bd2e40b;
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 6
// s needs to be of return_type, which is uint32_t
void seed(uint32_t s)
{
for (unsigned short i=0; i<4; ++i) {
_k[i] = 0;
_s[i] = 0;
}
_k[0] = s;
_o_counter = 0;
encrypt_counter();
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 7
template<class Seq>
void seed(Seq& q, typename sitmo_enable_if< has_generate_template<Seq>::value >::type* = 0 )
{
typename Seq::result_type w[8];
q.generate(&w[0], &w[8]);
for (unsigned short i=0; i<4; ++i) {
_k[i] = ( static_cast<uint64_t>(w[2*i]) << 32) | w[2*i+1];
_s[i] = 0;
}
_o_counter = 0;
encrypt_counter();
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 8
// Advances es state ei to ei+1 = TA(ei) and returns GA(ei).
uint32_t operator()()
{
// can we return a value from the current block?
if (_o_counter < 8) {
unsigned short _o_index = _o_counter >> 1;
_o_counter++;
if (_o_counter&1)
return _o[_o_index] & 0xFFFFFFFF;
else
return _o[_o_index] >> 32;
}
// generate a new block and return the first 32 bits
inc_counter();
encrypt_counter();
_o_counter = 1; // the next call
return _o[0] & 0xFFFFFFFF; // this call
}
// -------------------------------------------------
// misc
// -------------------------------------------------
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 9
// Advances es state ei to ei+z by any means equivalent to z
// consecutive calls e().
void discard(uint64_t z)
{
// check if we stay in the current block
if (z < 8 - _o_counter) {
_o_counter += static_cast<unsigned short>(z);
return;
}
// we will have to generate a new block...
z -= (8 - _o_counter); // discard the remainder of the current blok
_o_counter = z % 8; // set the pointer in the correct element in the new block
z -= _o_counter; // update z
z >>= 3; // the number of buffers is elements/8
++z; // and one more because we crossed the buffer line
inc_counter(z);
encrypt_counter();
}
// -------------------------------------------------
// IO
// -------------------------------------------------
template<class CharT, class Traits>
friend std::basic_ostream<CharT,Traits>&
operator<<(std::basic_ostream<CharT,Traits>& os, const prng_engine& s) {
for (unsigned short i=0; i<4; ++i)
os << s._k[i] << ' ' << s._s[i] << ' ' << s._o[i] << ' ';
os << s._o_counter;
return os;
}
template<class CharT, class Traits>
friend std::basic_istream<CharT,Traits>&
operator>>(std::basic_istream<CharT,Traits>& is, prng_engine& s) {
for (unsigned short i=0; i<4; ++i)
is >> s._k[i] >> s._s[i] >> s._o[i];
is >> s._o_counter;
return is;
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 10
// This operator is an equivalence relation. With Sx and Sy as the infinite
// sequences of values that would be generated by repeated future calls to
// x() and y(), respectively, returns true if Sx = Sy; else returns false.
bool operator==(const prng_engine& y)
{
if (_o_counter != y._o_counter) return false;
for (unsigned short i=0; i<4; ++i) {
if (_s[i] != y._s[i]) return false;
if (_k[i] != y._k[i]) return false;
if (_o[i] != y._o[i]) return false;
}
return true;
}
// req: 26.5.1.4 Random number engine requirements, p.908 table 117, row 11
bool operator!=(const prng_engine& y)
{
return !(*this == y);
}
// Extra function to set the key
void set_key(uint64_t k0=0, uint64_t k1=0, uint64_t k2=0, uint64_t k3=0)
{
_k[0] = k0; _k[1] = k1; _k[2] = k2; _k[3] = k3;
encrypt_counter();
}
// set the counter
void set_counter(uint64_t s0=0, uint64_t s1=0, uint64_t s2=0, uint64_t s3=0, unsigned short o_counter=0)
{
_s[0] = s0;
_s[1] = s1;
_s[2] = s2;
_s[3] = s3;
_o_counter = o_counter % 8;
encrypt_counter();
}
// versioning
uint32_t version()
{
return 5;
}
private:
void encrypt_counter()
{
uint64_t b[4];
uint64_t k[5];
for (unsigned short i=0; i<4; ++i) b[i] = _s[i];
for (unsigned short i=0; i<4; ++i) k[i] = _k[i];
k[4] = 0x1BD11BDAA9FC1A22 ^ k[0] ^ k[1] ^ k[2] ^ k[3];
MIXK(b[0], b[1], 14, b[2], b[3], 16, k[0], k[1], k[2], k[3]);
MIX2(b[0], b[3], 52, b[2], b[1], 57);
MIX2(b[0], b[1], 23, b[2], b[3], 40);
MIX2(b[0], b[3], 5, b[2], b[1], 37);
MIXK(b[0], b[1], 25, b[2], b[3], 33, k[1], k[2], k[3], k[4]+1);
MIX2(b[0], b[3], 46, b[2], b[1], 12);
MIX2(b[0], b[1], 58, b[2], b[3], 22);
MIX2(b[0], b[3], 32, b[2], b[1], 32);
MIXK(b[0], b[1], 14, b[2], b[3], 16, k[2], k[3], k[4], k[0]+2);
MIX2(b[0], b[3], 52, b[2], b[1], 57);
MIX2(b[0], b[1], 23, b[2], b[3], 40);
MIX2(b[0], b[3], 5, b[2], b[1], 37);
MIXK(b[0], b[1], 25, b[2], b[3], 33, k[3], k[4], k[0], k[1]+3);
MIX2(b[0], b[3], 46, b[2], b[1], 12);
MIX2(b[0], b[1], 58, b[2], b[3], 22);
MIX2(b[0], b[3], 32, b[2], b[1], 32);
MIXK(b[0], b[1], 14, b[2], b[3], 16, k[4], k[0], k[1], k[2]+4);
MIX2(b[0], b[3], 52, b[2], b[1], 57);
MIX2(b[0], b[1], 23, b[2], b[3], 40);
MIX2(b[0], b[3], 5, b[2], b[1], 37);
for (unsigned int i=0; i<4; ++i) _o[i] = b[i] + k[i];
_o[3] += 5;
}
void inc_counter()
{
++_s[0];
if (_s[0] != 0) return;
++_s[1];
if (_s[1] != 0) return;
++_s[2];
if (_s[2] != 0) return;
++_s[3];
}
void inc_counter(uint64_t z)
{
if (z > 0xFFFFFFFFFFFFFFFF - _s[0]) { // check if we will overflow the first 64 bit int
++_s[1];
if (_s[1] == 0) {
++_s[2];
if (_s[2] == 0) {
++_s[3];
}
}
}
_s[0] += z;
}
private:
uint64_t _k[4]; // key
uint64_t _s[4]; // state (counter)
uint64_t _o[4]; // cipher output 4 * 64 bit = 256 bit output
unsigned short _o_counter; // output chunk counter, the 256 random bits in _o
// are returned in eight 32 bit chunks
};
} // namespace sitmo
#undef MIXK
#undef MIX2
#endif

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
#include <algorithm> #include <algorithm>
namespace Grid { namespace Grid {

View File

@ -26,7 +26,7 @@ Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
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/Grid.h>
namespace Grid { namespace Grid {
} }

View File

@ -4,8 +4,9 @@ home=`pwd`
# library Make.inc # library Make.inc
cd $home/lib cd $home/lib
HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/Old/*' -not -path '*/Eigen/*'` HFILES=`find . -type f -name '*.h' -not -name '*Hdf5*' -not -path '*/gamma-gen/*' -not -path '*/Old/*' -not -path '*/Eigen/*'`
CCFILES=`find . -type f -name '*.cc' -not -name '*Communicator*.cc' -not -name '*Hdf5*'` HFILES="$HFILES"
CCFILES=`find . -type f -name '*.cc' -not -path '*/gamma-gen/*' -not -name '*Communicator*.cc' -not -name '*Hdf5*'`
echo HFILES=$HFILES > Make.inc echo HFILES=$HFILES > Make.inc
echo >> Make.inc echo >> Make.inc
echo CCFILES=$CCFILES >> Make.inc echo CCFILES=$CCFILES >> Make.inc

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
typedef DomainWallFermion<DomainWallVec5dImplR> DomainWallVecFermionR; typedef DomainWallFermion<DomainWallVec5dImplR> DomainWallVecFermionR;
@ -339,7 +339,7 @@ void TestMoo(This & Dw, That &sDw)
LatticeFermion ndiff(ngrid); LatticeFermion ndiff(ngrid);
LatticeFermion sdiff(sgrid); LatticeFermion sdiff(sgrid);
Gamma g5( Gamma::Gamma5 ); Gamma g5( Gamma::Algebra::Gamma5 );
std::vector<int> seeds({1,2,3,4,5,7,8}); std::vector<int> seeds({1,2,3,4,5,7,8});
GridParallelRNG RNG5(ngrid); GridParallelRNG RNG5(ngrid);

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int toint(const char* str){ int toint(const char* str){

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR; typedef WilsonFermion5D<DomainWallVec5dImplR> WilsonFermion5DR;

View File

@ -177,7 +177,7 @@ int main (int argc, char ** argv)
const int sdir=0; const int sdir=0;
RealD mass=0.01; RealD mass=0.01;
RealD M5 =1.0; RealD M5 =1.0;
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID); GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,&GRID);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,&GRID);
@ -218,12 +218,12 @@ int main (int argc, char ** argv)
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
// work out the predicted from Fourier // work out the predicted from Fourier
///////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT, Gamma::Algebra::GammaT,
Gamma::Gamma5 Gamma::Algebra::Gamma5
}; };
LatticeFermionD Kinetic(FGrid); Kinetic = zero; LatticeFermionD Kinetic(FGrid); Kinetic = zero;
LatticeComplexD kmu(FGrid); LatticeComplexD kmu(FGrid);
@ -311,7 +311,7 @@ int main (int argc, char ** argv)
std::cout << " Solving by FFT and Feynman rules" <<std::endl; std::cout << " Solving by FFT and Feynman rules" <<std::endl;
Ddwf.FreePropagator(src,ref,mass) ; Ddwf.FreePropagator(src,ref,mass) ;
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
LatticeFermionD src5(FGrid); src5=zero; LatticeFermionD src5(FGrid); src5=zero;
LatticeFermionD tmp5(FGrid); LatticeFermionD tmp5(FGrid);
@ -391,7 +391,7 @@ int main (int argc, char ** argv)
std::cout << " Solving by FFT and Feynman rules" <<std::endl; std::cout << " Solving by FFT and Feynman rules" <<std::endl;
Dov.FreePropagator(src,ref,mass) ; Dov.FreePropagator(src,ref,mass) ;
Gamma G5(Gamma::Gamma5); Gamma G5(Gamma::Algebra::Gamma5);
LatticeFermionD src5(FGrid); src5=zero; LatticeFermionD src5(FGrid); src5=zero;
LatticeFermionD tmp5(FGrid); LatticeFermionD tmp5(FGrid);

View File

@ -1,43 +1,218 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_gamma.cc Source file: ./tests/Test_gamma.cc
Copyright (C) 2015 Copyright (C) 2015-2017
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: Antonin Portelli <antonin.portelli@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/Grid.h> #include <Grid/Grid.h>
using namespace std; using namespace std;
using namespace Grid; using namespace Grid;
using namespace Grid::QCD; using namespace QCD;
//template<class vobj> class is_pod< iScalar<vobj> > static constexpr double tolerance = 1.0e-6;
//{ static std::array<SpinMatrix, Gamma::nGamma> testAlgebra;
//
//};
int main (int argc, char ** argv) void print(const SpinMatrix &g)
{
for(int i = 0; i < Ns; i++)
{
std::cout << GridLogMessage << "(";
for(int j=0;j<Ns;j++){
if ( abs(g()(i,j)()) == 0 ) {
std::cout<< " 0";
} else if ( abs(g()(i,j)() - Complex(0,1)) == 0){
std::cout<< " i";
} else if ( abs(g()(i,j)() + Complex(0,1)) == 0){
std::cout<< "-i";
} else if ( abs(g()(i,j)() - Complex(1,0)) == 0){
std::cout<< " 1";
} else if ( abs(g()(i,j)() + Complex(1,0)) == 0){
std::cout<< "-1";
}
std::cout<<((j == Ns-1) ? ")" : "," );
}
std::cout << std::endl;
}
std::cout << GridLogMessage << std::endl;
}
void createTestAlgebra(void)
{
std::array<SpinMatrix, 4> testg;
SpinMatrix testg5;
const Complex I(0., 1.), mI(0., -1.);
testg[0] = zero;
testg[0]()(0, 3) = I;
testg[0]()(1, 2) = I;
testg[0]()(2, 1) = mI;
testg[0]()(3, 0) = mI;
std::cout << GridLogMessage << "test GammaX= " << std::endl;
print(testg[0]);
testg[1] = zero;
testg[1]()(0, 3) = -1.;
testg[1]()(1, 2) = 1.;
testg[1]()(2, 1) = 1.;
testg[1]()(3, 0) = -1.;
std::cout << GridLogMessage << "test GammaY= " << std::endl;
print(testg[1]);
testg[2] = zero;
testg[2]()(0, 2) = I;
testg[2]()(1, 3) = mI;
testg[2]()(2, 0) = mI;
testg[2]()(3, 1) = I;
std::cout << GridLogMessage << "test GammaZ= " << std::endl;
print(testg[2]);
testg[3] = zero;
testg[3]()(0, 2) = 1.;
testg[3]()(1, 3) = 1.;
testg[3]()(2, 0) = 1.;
testg[3]()(3, 1) = 1.;
std::cout << GridLogMessage << "test GammaT= " << std::endl;
print(testg[3]);
testg5 = testg[0]*testg[1]*testg[2]*testg[3];
#define DEFINE_TEST_G(g, exp)\
testAlgebra[Gamma::Algebra::g] = exp;\
testAlgebra[Gamma::Algebra::Minus##g] = -exp;\
DEFINE_TEST_G(Identity , 1.);
DEFINE_TEST_G(Gamma5 , testg5);
DEFINE_TEST_G(GammaX , testg[0]);
DEFINE_TEST_G(GammaY , testg[1]);
DEFINE_TEST_G(GammaZ , testg[2]);
DEFINE_TEST_G(GammaT , testg[3]);
DEFINE_TEST_G(GammaXGamma5, testg[0]*testg5);
DEFINE_TEST_G(GammaYGamma5, testg[1]*testg5);
DEFINE_TEST_G(GammaZGamma5, testg[2]*testg5);
DEFINE_TEST_G(GammaTGamma5, testg[3]*testg5);
DEFINE_TEST_G(SigmaXY , .5*(testg[0]*testg[1] - testg[1]*testg[0]));
DEFINE_TEST_G(SigmaXZ , .5*(testg[0]*testg[2] - testg[2]*testg[0]));
DEFINE_TEST_G(SigmaXT , .5*(testg[0]*testg[3] - testg[3]*testg[0]));
DEFINE_TEST_G(SigmaYZ , .5*(testg[1]*testg[2] - testg[2]*testg[1]));
DEFINE_TEST_G(SigmaYT , .5*(testg[1]*testg[3] - testg[3]*testg[1]));
DEFINE_TEST_G(SigmaZT , .5*(testg[2]*testg[3] - testg[3]*testg[2]));
#undef DEFINE_TEST_G
}
template <typename Expr>
void test(const Expr &a, const Expr &b)
{
if (norm2(a - b) < tolerance)
{
std::cout << "[OK] ";
}
else
{
std::cout << "[fail]" << std::endl;
std::cout << GridLogError << "a= " << a << std::endl;
std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
std::cout << GridLogError << "b= " << b << std::endl;
exit(EXIT_FAILURE);
}
}
void checkMat(const Gamma::Algebra a, GridSerialRNG &rng)
{
SpinVector v;
SpinMatrix m, &testg = testAlgebra[a];
Gamma g(a);
bool pass = true;
random(rng, v);
random(rng, m);
std::cout << GridLogMessage << "Checking " << Gamma::name[a] << ": ";
std::cout << "vecmul ";
test(g*v, testg*v);
std::cout << "matlmul ";
test(g*m, testg*m);
std::cout << "matrmul ";
test(m*g, m*testg);
std::cout << std::endl;
}
void checkProd(const Gamma::Algebra a, const Gamma::Algebra b)
{
SpinMatrix gm, testg = testAlgebra[a]*testAlgebra[b];
Gamma g = Gamma(a)*Gamma(b);
bool pass = true;
std::cout << GridLogMessage << "Checking " << Gamma::name[a] << " * "
<< Gamma::name[b] << ": ";
gm = 1.0;
gm = g*gm;
test(gm, testg);
std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
}
void checkAdj(const Gamma::Algebra a)
{
SpinMatrix gm, testg = adj(testAlgebra[a]);
Gamma g(adj(Gamma(a)));
bool pass = true;
std::cout << GridLogMessage << "Checking adj(" << Gamma::name[a] << "): ";
gm = 1.0;
gm = g*gm;
test(gm, testg);
std::cout << "(= " << Gamma::name[g.g] << ")" << std::endl;
}
void checkProject(GridSerialRNG &rng)
{
SpinVector rv, recon, full;
HalfSpinVector hsp, hsm;
random(rng, rv);
#define CHECK_PROJ(dir, gamma)\
std::cout << GridLogMessage << "Checking " << #dir << " projector: ";\
spProj##dir(hsm,rv);\
spRecon##dir(recon,hsm);\
test(recon, rv + Gamma(Gamma::Algebra::gamma)*rv);\
std::cout << std::endl;
CHECK_PROJ(Xp, GammaX);
CHECK_PROJ(Yp, GammaY);
CHECK_PROJ(Zp, GammaZ);
CHECK_PROJ(Tp, GammaT);
CHECK_PROJ(5p, Gamma5);
CHECK_PROJ(Xm, MinusGammaX);
CHECK_PROJ(Ym, MinusGammaY);
CHECK_PROJ(Zm, MinusGammaZ);
CHECK_PROJ(Tm, MinusGammaT);
CHECK_PROJ(5m, MinusGamma5);
#undef CHECK_PROJ
}
int main(int argc, char *argv[])
{ {
Grid_init(&argc,&argv); Grid_init(&argc,&argv);
@ -46,178 +221,36 @@ int main (int argc, char ** argv)
std::vector<int> mpi_layout = GridDefaultMpi(); std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout); GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridParallelRNG pRNG(&Grid);
pRNG.SeedRandomDevice();
GridSerialRNG sRNG; GridSerialRNG sRNG;
sRNG.SeedRandomDevice(); sRNG.SeedRandomDevice();
SpinMatrix ident; ident=zero; std::cout << GridLogMessage << "======== Test algebra" << std::endl;
SpinMatrix rnd ; random(sRNG,rnd); createTestAlgebra();
std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
SpinMatrix ll; ll=zero; for (int i = 0; i < Gamma::nGamma; ++i)
SpinMatrix rr; rr=zero; {
SpinMatrix result; checkMat(i, sRNG);
SpinVector lv; random(sRNG,lv);
SpinVector rv; random(sRNG,rv);
// std::cout<<GridLogMessage << " Is pod " << std::is_pod<SpinVector>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod double " << std::is_pod<double>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod ComplexF " << std::is_pod<ComplexF>::value << std::endl;
// std::cout<<GridLogMessage << " Is triv double " << std::has_trivial_default_constructor<double>::value << std::endl;
// std::cout<<GridLogMessage << " Is triv ComplexF " << std::has_trivial_default_constructor<ComplexF>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<double> " << std::is_pod<iScalar<double> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<ComplexF> " << std::is_pod<iScalar<ComplexF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vComplexF> " << std::is_pod<iScalar<vComplexF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vComplexD> " << std::is_pod<iScalar<vComplexD> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vRealF> " << std::is_pod<iScalar<vRealF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vRealD> " << std::is_pod<iScalar<vRealD> >::value << std::endl;
// std::cout<<GridLogMessage << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
// std::cout<<GridLogMessage << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value << std::endl;
for(int a=0;a<Ns;a++){
ident()(a,a) = ComplexF(1.0);
} }
std::cout << GridLogMessage << std::endl;
const Gamma::GammaMatrix *g = Gamma::GammaMatrices; std::cout << GridLogMessage << "======== Algebra multiplication table check" << std::endl;
const char **list = Gamma::GammaMatrixNames; for (int i = 0; i < Gamma::nGamma; ++i)
for (int j = 0; j < Gamma::nGamma; ++j)
result =ll*Gamma(g[0])*rr; {
result =ll*Gamma(g[0]); checkProd(i, j);
rv = Gamma(g[0])*lv;
for(int mu=0;mu<12;mu++){
result = Gamma(g[mu])* ident;
for(int i=0;i<Ns;i++){
if(i==0) std::cout<<GridLogMessage << list[mu];
else std::cout<<GridLogMessage << list[12];
std::cout<<"(";
for(int j=0;j<Ns;j++){
if ( abs(result()(i,j)())==0 ) {
std::cout<< " 0";
} else if ( abs(result()(i,j)() - Complex(0,1))==0){
std::cout<< " i";
} else if ( abs(result()(i,j)() + Complex(0,1))==0){
std::cout<< "-i";
} else if ( abs(result()(i,j)() - Complex(1,0))==0){
std::cout<< " 1";
} else if ( abs(result()(i,j)() + Complex(1,0))==0){
std::cout<< "-1";
} }
std::cout<<((j==Ns-1) ? ")" : "," ); std::cout << GridLogMessage << std::endl;
std::cout << GridLogMessage << "======== Adjoints check" << std::endl;
for (int i = 0; i < Gamma::nGamma; ++i)
{
checkAdj(i);
} }
std::cout << std::endl; std::cout << GridLogMessage << std::endl;
} std::cout << GridLogMessage << "======== Spin projectors check" << std::endl;
checkProject(sRNG);
std::cout << std::endl; std::cout << GridLogMessage << std::endl;
}
std::cout << "Testing Gamma^2 - 1 = 0"<<std::endl;
for(int mu=0;mu<6;mu++){
result = Gamma(g[mu])* ident * Gamma(g[mu]);
result = result - ident;
RealD mag = norm2(result);
std::cout << list[mu]<<" " << mag<<std::endl;
}
std::cout << "Testing (MinusGamma + G )M = 0"<<std::endl;
for(int mu=0;mu<6;mu++){
result = rnd * Gamma(g[mu]);
result = result + rnd * Gamma(g[mu+6]);
RealD mag = norm2(result);
std::cout << list[mu]<<" " << mag<<std::endl;
}
std::cout << "Testing M(MinusGamma + G ) = 0"<<std::endl;
for(int mu=0;mu<6;mu++){
result = Gamma(g[mu]) *rnd;
result = result + Gamma(g[mu+6])*rnd;
RealD mag = norm2(result);
std::cout << list[mu]<<" " << mag<<std::endl;
}
// Testing spins and reconstructs
SpinVector recon; random(sRNG,rv);
SpinVector full;
HalfSpinVector hsp,hsm;
// Xp
double mag;
spProjXp(hsm,rv);
spReconXp(recon,hsm);
full = rv + Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Xp "<< mag<<std::endl;
// Xm
spProjXm(hsm,rv);
spReconXm(recon,hsm);
full = rv - Gamma(Gamma::GammaX) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Xm "<< mag<<std::endl;
// Yp
spProjYp(hsm,rv);
spReconYp(recon,hsm);
full = rv + Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Yp "<< mag<<std::endl;
// Ym
spProjYm(hsm,rv);
spReconYm(recon,hsm);
full = rv - Gamma(Gamma::GammaY) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Ym "<< mag<<std::endl;
// Zp
spProjZp(hsm,rv);
spReconZp(recon,hsm);
full = rv + Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Zp "<< mag<<std::endl;
// Zm
spProjZm(hsm,rv);
spReconZm(recon,hsm);
full = rv - Gamma(Gamma::GammaZ) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Zm "<< mag<<std::endl;
// Tp
spProjTp(hsm,rv);
spReconTp(recon,hsm);
full = rv + Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Tp "<< mag<<std::endl;
// Tm
spProjTm(hsm,rv);
spReconTm(recon,hsm);
full = rv - Gamma(Gamma::GammaT) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "Tm "<< mag<<std::endl;
// 5p
spProj5p(hsm,rv);
spRecon5p(recon,hsm);
full = rv + Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "5p "<< mag<<std::endl;
// 5m
spProj5m(hsm,rv);
spRecon5m(recon,hsm);
full = rv - Gamma(Gamma::Gamma5) *rv;
mag = TensorRemove(norm2(full-recon));
std::cout << "5m "<< mag<<std::endl;
Grid_finalize(); Grid_finalize();
return EXIT_SUCCESS;
} }

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
template<class What> template<class What>

View File

@ -38,11 +38,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };

View File

@ -131,6 +131,7 @@ int main(int argc, char *argv[])
mesPar.q1 = qName[i]; mesPar.q1 = qName[i];
mesPar.q2 = qName[j]; mesPar.q2 = qName[j];
mesPar.gammas = "all"; mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_Z2_" application.createModule<MContraction::Meson>("meson_Z2_"
+ std::to_string(t) + std::to_string(t)
+ "_" + "_"
@ -149,6 +150,7 @@ int main(int argc, char *argv[])
mesPar.q1 = qName[i]; mesPar.q1 = qName[i];
mesPar.q2 = seqName[j][mu]; mesPar.q2 = seqName[j][mu];
mesPar.gammas = "all"; mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("3pt_Z2_" application.createModule<MContraction::Meson>("3pt_Z2_"
+ std::to_string(t) + std::to_string(t)
+ "_" + "_"

View File

@ -97,6 +97,7 @@ int main(int argc, char *argv[])
mesPar.q1 = "Qpt_" + flavour[i]; mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j]; mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "all"; mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_pt_" application.createModule<MContraction::Meson>("meson_pt_"
+ flavour[i] + flavour[j], + flavour[i] + flavour[j],
mesPar); mesPar);
@ -104,6 +105,7 @@ int main(int argc, char *argv[])
mesPar.q1 = "QZ2_" + flavour[i]; mesPar.q1 = "QZ2_" + flavour[i];
mesPar.q2 = "QZ2_" + flavour[j]; mesPar.q2 = "QZ2_" + flavour[j];
mesPar.gammas = "all"; mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_Z2_" application.createModule<MContraction::Meson>("meson_Z2_"
+ flavour[i] + flavour[j], + flavour[i] + flavour[j],
mesPar); mesPar);

View File

@ -25,7 +25,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
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/Grid.h>
double calc_grid_p (Grid::QCD::LatticeGaugeField & lat); double calc_grid_p (Grid::QCD::LatticeGaugeField & lat);
double calc_chroma_p (Grid::QCD::LatticeGaugeField & lat); double calc_chroma_p (Grid::QCD::LatticeGaugeField & lat);

View File

@ -26,7 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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/Grid.h>
int Ls=8; int Ls=8;
double M5=1.6; double M5=1.6;

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };

View File

@ -37,8 +37,8 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu[] = {Gamma::GammaX, Gamma::GammaY, Gamma::GammaZ, Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ,
Gamma::GammaT}; Gamma::Algebra::GammaT};
int main(int argc, char** argv) { int main(int argc, char** argv) {
Grid_init(&argc, &argv); Grid_init(&argc, &argv);

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -38,11 +38,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -504,7 +504,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds); GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
Gamma g5(Gamma::Gamma5); Gamma g5(Gamma::Algebra::Gamma5);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src; LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid); result=zero; LatticeFermion result(FGrid); result=zero;

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -37,11 +37,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)

View File

@ -36,11 +36,11 @@ struct scal {
d internal; d internal;
}; };
Gamma::GammaMatrix Gmu [] = { Gamma::Algebra Gmu [] = {
Gamma::GammaX, Gamma::Algebra::GammaX,
Gamma::GammaY, Gamma::Algebra::GammaY,
Gamma::GammaZ, Gamma::Algebra::GammaZ,
Gamma::GammaT Gamma::Algebra::GammaT
}; };
int main (int argc, char ** argv) int main (int argc, char ** argv)