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

Merge branch 'develop' into feature/hmc_generalise

Browse Source
This commit is contained in:
Guido Cossu
2017-02-09 16:12:00 +00:00
parent c67f41887b c94133af49
commit e0571c872b
85 changed files with 5136 additions and 1145 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
lib/Init.cc
View File

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

2
lib/Log.cc
View File

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

4
lib/PerfCount.cc
View File

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

6
lib/Stat.cc
View File

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

2
lib/algorithms/approx/MultiShiftFunction.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid {
double MultiShiftFunction::approx(double x)

2
lib/algorithms/approx/Remez.cc
View File

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

9
lib/algorithms/iterative/ConjugateGradient.h
View File

@@ -45,6 +45,8 @@ class ConjugateGradient : public OperatorFunction<Field> {
// Defaults true.
RealD Tolerance;
Integer MaxIterations;
Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true)
: Tolerance(tol),
MaxIterations(maxit),
@@ -158,13 +160,16 @@ class ConjugateGradient : public OperatorFunction<Field> {
std::cout << std::endl;
if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
IterationsToComplete = k;
return;
}
}
std::cout << GridLogMessage << "ConjugateGradient did NOT converge"
<< std::endl;
if (ErrorOnNoConverge) exit(1);
if (ErrorOnNoConverge) assert(0);
IterationsToComplete = k;
}
};
}

20
lib/algorithms/iterative/ConjugateGradientMixedPrec.h
View File

@@ -35,6 +35,7 @@ namespace Grid {
class MixedPrecisionConjugateGradient : public LinearFunction<FieldD> {
public:
RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations;
Integer MaxOuterIterations;
GridBase* SinglePrecGrid; //Grid for single-precision fields
@@ -42,12 +43,16 @@ namespace Grid {
LinearOperatorBase<FieldF> &Linop_f;
LinearOperatorBase<FieldD> &Linop_d;
Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
//Option to speed up *inner single precision* solves using a LinearFunction that produces a guess
LinearFunction<FieldF> *guesser;
MixedPrecisionConjugateGradient(RealD tol, Integer maxinnerit, Integer maxouterit, GridBase* _sp_grid, LinearOperatorBase<FieldF> &_Linop_f, LinearOperatorBase<FieldD> &_Linop_d) :
Linop_f(_Linop_f), Linop_d(_Linop_d),
Tolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid),
Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid),
OuterLoopNormMult(100.), guesser(NULL){ };
void useGuesser(LinearFunction<FieldF> &g){
@@ -55,6 +60,8 @@ namespace Grid {
}
void operator() (const FieldD &src_d_in, FieldD &sol_d){
TotalInnerIterations = 0;
GridStopWatch TotalTimer;
TotalTimer.Start();
@@ -74,7 +81,7 @@ namespace Grid {
FieldD src_d(DoublePrecGrid);
src_d = src_d_in; //source for next inner iteration, computed from residual during operation
RealD inner_tol = Tolerance;
RealD inner_tol = InnerTolerance;
FieldF src_f(SinglePrecGrid);
src_f.checkerboard = cb;
@@ -89,7 +96,9 @@ namespace Grid {
GridStopWatch PrecChangeTimer;
for(Integer outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){
Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count
for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++){
//Compute double precision rsd and also new RHS vector.
Linop_d.HermOp(sol_d, tmp_d);
RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector
@@ -117,6 +126,7 @@ namespace Grid {
InnerCGtimer.Start();
CG_f(Linop_f, src_f, sol_f);
InnerCGtimer.Stop();
TotalInnerIterations += CG_f.IterationsToComplete;
//Convert sol back to double and add to double prec solution
PrecChangeTimer.Start();
@@ -131,9 +141,11 @@ namespace Grid {
ConjugateGradient<FieldD> CG_d(Tolerance, MaxInnerIterations);
CG_d(Linop_d, src_d_in, sol_d);
TotalFinalStepIterations = CG_d.IterationsToComplete;
TotalTimer.Stop();
std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Total " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl;
std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl;
std::cout<<GridLogMessage<<"MixedPrecisionConjugateGradient: Total time " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl;
}
};

2
lib/algorithms/iterative/Matrix.h
View File

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

2
lib/communicator/Communicator_base.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
namespace Grid {
///////////////////////////////////////////////////////////////

2
lib/communicator/Communicator_mpi.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
#include <mpi.h>
namespace Grid {

2
lib/communicator/Communicator_mpi3.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
#include <mpi.h>
namespace Grid {

2
lib/communicator/Communicator_none.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
namespace Grid {
///////////////////////////////////////////////////////////////////////////////////////////////////

2
lib/communicator/Communicator_shmem.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
#include <mpp/shmem.h>
namespace Grid {

7
lib/lattice/Lattice_rng.h
View File

@@ -30,6 +30,7 @@
#define GRID_LATTICE_RNG_H
#include <random>
#include <Grid/sitmo_rng/sitmo_prng_engine.hpp>
namespace Grid {
@@ -132,10 +133,14 @@ namespace Grid {
typedef uint64_t RngStateType;
typedef std::ranlux48 RngEngine;
static const int RngStateCount = 15;
#else
#elif RNG_MT19937
typedef std::mt19937 RngEngine;
typedef uint32_t RngStateType;
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
std::vector<RngEngine> _generators;
std::vector<std::uniform_real_distribution<RealD>> _uniform;

1
lib/lattice/Lattice_transfer.h
View File

@@ -386,6 +386,7 @@ void InsertSlice(Lattice<vobj> &lowDim,Lattice<vobj> & higherDim,int slice, int
}
// the above should guarantee that the operations are local
// Guido: check the threading here
//PARALLEL_FOR_LOOP
for(int idx=0;idx<lg->lSites();idx++){
std::vector<int> lcoor(nl);

2
lib/pugixml/pugixml.cc
View File

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

2
lib/qcd/action/fermion/CayleyFermion5D.cc
View File

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

2
lib/qcd/action/fermion/CayleyFermion5Dcache.cc
View File

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

2
lib/qcd/action/fermion/CayleyFermion5Ddense.cc
View File

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

2
lib/qcd/action/fermion/CayleyFermion5Dssp.cc
View File

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

2
lib/qcd/action/fermion/CayleyFermion5Dvec.cc
View File

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

2
lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid {
namespace QCD {

2
lib/qcd/action/fermion/PartialFractionFermion5D.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid {
namespace QCD {

12
lib/qcd/action/fermion/WilsonFermion.cc
View File

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

24
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@@ -30,8 +30,8 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <PerfCount.h>
#include <Grid/Grid.h>
#include <Grid/PerfCount.h>
namespace Grid {
namespace QCD {
@@ -464,11 +464,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
typedef iSinglet<ScalComplex> Tcomplex;
typedef Lattice<iSinglet<vector_type> > LatComplex;
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
Gamma::Algebra Gmu [] = {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT
};
std::vector<int> latt_size = _grid->_fdimensions;
@@ -535,11 +535,11 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
template<class Impl>
void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass)
{
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
Gamma::Algebra Gmu [] = {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT
};
GridBase *_grid = _FourDimGrid;

2
lib/qcd/action/fermion/WilsonKernels.cc
View File

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

2
lib/qcd/action/fermion/WilsonKernelsAsm.cc
View File

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

2
lib/qcd/action/fermion/WilsonKernelsHand.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
#define REGISTER

2
lib/qcd/action/fermion/WilsonTMFermion.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid {
namespace QCD {

2
lib/qcd/action/fermion/g5HermitianLinop.h
View File

@@ -80,7 +80,7 @@ class Gamma5HermitianLinearOperator : public LinearOperatorBase<Field> {
Matrix &_Mat;
Gamma g5;
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){
HermOp(in,out);
}

2
lib/qcd/hmc/HMC.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid{
namespace QCD{

95
lib/qcd/spin/Dirac.cc
View File

@@ -1,95 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/spin/Dirac.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
namespace Grid {
namespace QCD {
Gamma::GammaMatrix Gamma::GammaMatrices [] = {
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 ",
" "
};
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);
// }
}
}

796
lib/qcd/spin/Dirac.h
View File

@@ -1,628 +1,232 @@
/*************************************************************************************
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/spin/Dirac.h
Source file: lib/qcd/spin/Dirac.h
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 <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef GRID_QCD_DIRAC_H
#define GRID_QCD_DIRAC_H
namespace Grid{
// Gamma matrices using the code generated by the Mathematica notebook
// gamma-gen/gamma-gen.nb in Gamma.cc & Gamma.h
////////////////////////////////////////////////////////////////////////////////
#include <Grid/qcd/spin/Gamma.h>
namespace Grid {
// Dirac algebra adjoint operator (not in QCD:: to overload other adj)
inline QCD::Gamma adj(const QCD::Gamma &g)
{
return QCD::Gamma (QCD::Gamma::adj[g.g]);
}
namespace QCD {
// Dirac algebra mutliplication operator
inline Gamma operator*(const Gamma &g1, const Gamma &g2)
{
return Gamma (Gamma::mul[g1.g][g2.g]);
}
class Gamma {
// general left multiply
template<class vtype>
inline auto operator*(const Gamma &G, const iScalar<vtype> &arg)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type
{
iScalar<vtype> ret;
ret._internal=G*arg._internal;
return ret;
}
public:
template<class vtype,int N>
inline auto operator*(const Gamma &G, const iVector<vtype, N> &arg)
->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i]=G*arg._internal[i];
}
return ret;
}
const int Ns=4;
enum GammaMatrix {
Identity,
GammaX,
GammaY,
GammaZ,
GammaT,
Gamma5,
MinusIdentity,
MinusGammaX,
MinusGammaY,
MinusGammaZ,
MinusGammaT,
MinusGamma5
// GammaXGamma5, // Rest are composite (willing to take hit for two calls sequentially)
// GammaYGamma5, // as they are less commonly used.
// GammaZGamma5,
// GammaTGamma5,
// SigmaXY,
// SigmaXZ,
// SigmaYZ,
// SigmaXT,
// SigmaYT,
// SigmaZT,
// MinusGammaXGamma5, easiest to form by composition
// MinusGammaYGamma5, as performance is not critical for these
// MinusGammaZGamma5,
// MinusGammaTGamma5,
// MinusSigmaXY,
// MinusSigmaXZ,
// MinusSigmaYZ,
// MinusSigmaXT,
// MinusSigmaYT,
// MinusSigmaZT
};
static GammaMatrix GammaMatrices[];
static const char *GammaMatrixNames[];
template<class vtype, int N>
inline auto operator*(const Gamma &G, const iMatrix<vtype, N> &arg)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=G*arg._internal[i][j];
}}
return ret;
}
Gamma (GammaMatrix g) { _g=g; }
// general right multiply
template<class vtype>
inline auto operator*(const iScalar<vtype> &arg, const Gamma &G)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type
{
iScalar<vtype> ret;
ret._internal=arg._internal*G;
return ret;
}
GammaMatrix _g;
template<class vtype, int N>
inline auto operator * (const iMatrix<vtype, N> &arg, const Gamma &G)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=arg._internal[i][j]*G;
}}
return ret;
}
};
// Gamma-left matrices gL_mu = g_mu*(1 - g5)
////////////////////////////////////////////////////////////////////////////////
class GammaL
{
public:
typedef Gamma::Algebra Algebra;
Gamma gamma;
public:
GammaL(const Algebra initg): gamma(initg) {}
GammaL(const Gamma initg): gamma(initg) {}
};
// vector multiply
template<class vtype>
inline auto operator*(const GammaL &gl, const iVector<vtype, Ns> &arg)
->typename std::enable_if<matchGridTensorIndex<iVector<vtype, Ns>, SpinorIndex>::value, iVector<vtype, Ns>>::type
{
iVector<vtype, Ns> buf;
// Make gamma products (Chroma convention)
SpinMatrix makeGammaProd(const unsigned int i);
/* Gx
* 0 0 0 i
* 0 0 i 0
* 0 -i 0 0
* -i 0 0 0
*/
// right multiplication makes sense for matrix args, not for vector since there is
// no concept of row versus columnar indices
template<class vtype> inline void rmultMinusGammaX(iMatrix<vtype,Ns> &ret,const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = timesI(rhs(i,3));
ret(i,1) = timesI(rhs(i,2));
ret(i,2) = timesMinusI(rhs(i,1));
ret(i,3) = timesMinusI(rhs(i,0));
}
};
template<class vtype> inline void rmultGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = timesMinusI(rhs(i,3));
ret(i,1) = timesMinusI(rhs(i,2));
ret(i,2) = timesI(rhs(i,1));
ret(i,3) = timesI(rhs(i,0));
}
};
template<class vtype> inline void multGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = timesI(rhs(3,i));
ret(1,i) = timesI(rhs(2,i));
ret(2,i) = timesMinusI(rhs(1,i));
ret(3,i) = timesMinusI(rhs(0,i));
}
};
template<class vtype> inline void multMinusGammaX(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = timesMinusI(rhs(3,i));
ret(1,i) = timesMinusI(rhs(2,i));
ret(2,i) = timesI(rhs(1,i));
ret(3,i) = timesI(rhs(0,i));
}
};
buf(0) = 0.;
buf(1) = 0.;
buf(2) = 2.*arg(2);
buf(3) = 2.*arg(3);
return gl.gamma*buf;
};
template<class vtype> inline void multGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret._internal[0] = timesI(rhs._internal[3]);
ret._internal[1] = timesI(rhs._internal[2]);
ret._internal[2] = timesMinusI(rhs._internal[1]);
ret._internal[3] = timesMinusI(rhs._internal[0]);
};
template<class vtype> inline void multMinusGammaX(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = timesMinusI(rhs(3));
ret(1) = timesMinusI(rhs(2));
ret(2) = timesI(rhs(1));
ret(3) = timesI(rhs(0));
};
// matrix left multiply
template<class vtype>
inline auto operator*(const GammaL &gl, const iMatrix<vtype, Ns> &arg)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
{
iMatrix<vtype, Ns> buf;
for(unsigned int i = 0; i < Ns; ++i)
{
buf(0, i) = 0.;
buf(1, i) = 0.;
buf(2, i) = 2.*arg(2, i);
buf(3, i) = 2.*arg(3, i);
}
return gl.gamma*buf;
};
// matrix right multiply
template<class vtype>
inline auto operator*(const iMatrix<vtype, Ns> &arg, const GammaL &gl)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype, Ns>, SpinorIndex>::value, iMatrix<vtype, Ns>>::type
{
iMatrix<vtype, Ns> buf;
buf = arg*gl.gamma;
for(unsigned int i = 0; i < Ns; ++i)
{
buf(i, 0) = 0.;
buf(i, 1) = 0.;
buf(i, 2) = 2.*buf(i, 2);
buf(i, 3) = 2.*buf(i, 3);
}
return buf;
};
/*Gy
* 0 0 0 -1 [0] -+ [3]
* 0 0 1 0 [1] +- [2]
* 0 1 0 0
* -1 0 0 0
*/
template<class vtype> inline void rmultGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = -rhs(i,3);
ret(i,1) = rhs(i,2);
ret(i,2) = rhs(i,1);
ret(i,3) = -rhs(i,0);
}
};
template<class vtype> inline void rmultMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = rhs(i,3);
ret(i,1) = -rhs(i,2);
ret(i,2) = -rhs(i,1);
ret(i,3) = rhs(i,0);
}
};
template<class vtype> inline void multGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = -rhs(3,i);
ret(1,i) = rhs(2,i);
ret(2,i) = rhs(1,i);
ret(3,i) = -rhs(0,i);
}
};
template<class vtype> inline void multMinusGammaY(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = rhs(3,i);
ret(1,i) = -rhs(2,i);
ret(2,i) = -rhs(1,i);
ret(3,i) = rhs(0,i);
}
};
template<class vtype> inline void multGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = -rhs(3);
ret(1) = rhs(2);
ret(2) = rhs(1);
ret(3) = -rhs(0);
};
template<class vtype> inline void multMinusGammaY(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = rhs(3);
ret(1) = -rhs(2);
ret(2) = -rhs(1);
ret(3) = rhs(0);
};
/*Gz
* 0 0 i 0 [0]+-i[2]
* 0 0 0 -i [1]-+i[3]
* -i 0 0 0
* 0 i 0 0
*/
template<class vtype> inline void rmultGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = timesMinusI(rhs(i,2));
ret(i,1) = timesI(rhs(i,3));
ret(i,2) = timesI(rhs(i,0));
ret(i,3) = timesMinusI(rhs(i,1));
}
};
template<class vtype> inline void rmultMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = timesI(rhs(i,2));
ret(i,1) = timesMinusI(rhs(i,3));
ret(i,2) = timesMinusI(rhs(i,0));
ret(i,3) = timesI(rhs(i,1));
}
};
template<class vtype> inline void multGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = timesI(rhs(2,i));
ret(1,i) =timesMinusI(rhs(3,i));
ret(2,i) =timesMinusI(rhs(0,i));
ret(3,i) = timesI(rhs(1,i));
}
};
template<class vtype> inline void multMinusGammaZ(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = timesMinusI(rhs(2,i));
ret(1,i) = timesI(rhs(3,i));
ret(2,i) = timesI(rhs(0,i));
ret(3,i) = timesMinusI(rhs(1,i));
}
};
template<class vtype> inline void multGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = timesI(rhs(2));
ret(1) =timesMinusI(rhs(3));
ret(2) =timesMinusI(rhs(0));
ret(3) = timesI(rhs(1));
};
template<class vtype> inline void multMinusGammaZ(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = timesMinusI(rhs(2));
ret(1) = timesI(rhs(3));
ret(2) = timesI(rhs(0));
ret(3) = timesMinusI(rhs(1));
};
/*Gt
* 0 0 1 0 [0]+-[2]
* 0 0 0 1 [1]+-[3]
* 1 0 0 0
* 0 1 0 0
*/
template<class vtype> inline void rmultGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = rhs(i,2);
ret(i,1) = rhs(i,3);
ret(i,2) = rhs(i,0);
ret(i,3) = rhs(i,1);
}
};
template<class vtype> inline void rmultMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) =- rhs(i,2);
ret(i,1) =- rhs(i,3);
ret(i,2) =- rhs(i,0);
ret(i,3) =- rhs(i,1);
}
};
template<class vtype> inline void multGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = rhs(2,i);
ret(1,i) = rhs(3,i);
ret(2,i) = rhs(0,i);
ret(3,i) = rhs(1,i);
}
};
template<class vtype> inline void multMinusGammaT(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) =-rhs(2,i);
ret(1,i) =-rhs(3,i);
ret(2,i) =-rhs(0,i);
ret(3,i) =-rhs(1,i);
}
};
template<class vtype> inline void multGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = rhs(2);
ret(1) = rhs(3);
ret(2) = rhs(0);
ret(3) = rhs(1);
};
template<class vtype> inline void multMinusGammaT(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) =-rhs(2);
ret(1) =-rhs(3);
ret(2) =-rhs(0);
ret(3) =-rhs(1);
};
/*G5
* 1 0 0 0 [0]+-[2]
* 0 1 0 0 [1]+-[3]
* 0 0 -1 0
* 0 0 0 -1
*/
template<class vtype> inline void rmultGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) = rhs(i,0);
ret(i,1) = rhs(i,1);
ret(i,2) =-rhs(i,2);
ret(i,3) =-rhs(i,3);
}
};
template<class vtype> inline void rmultMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(i,0) =-rhs(i,0);
ret(i,1) =-rhs(i,1);
ret(i,2) = rhs(i,2);
ret(i,3) = rhs(i,3);
}
};
//general left multiply
template<class vtype>
inline auto operator*(const GammaL &gl, const iScalar<vtype> &arg)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type
{
iScalar<vtype> ret;
ret._internal=gl*arg._internal;
return ret;
}
template<class vtype> inline void multGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) = rhs(0,i);
ret(1,i) = rhs(1,i);
ret(2,i) =-rhs(2,i);
ret(3,i) =-rhs(3,i);
}
};
template<class vtype> inline void multMinusGamma5(iMatrix<vtype,Ns> &ret, const iMatrix<vtype,Ns> &rhs){
for(int i=0;i<Ns;i++){
ret(0,i) =-rhs(0,i);
ret(1,i) =-rhs(1,i);
ret(2,i) = rhs(2,i);
ret(3,i) = rhs(3,i);
}
};
template<class vtype,int N>
inline auto operator*(const GammaL &gl, const iVector<vtype, N> &arg)
->typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N>>::type
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i]=gl*arg._internal[i];
}
return ret;
}
template<class vtype> inline void multGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) = rhs(0);
ret(1) = rhs(1);
ret(2) =-rhs(2);
ret(3) =-rhs(3);
};
template<class vtype> inline void multMinusGamma5(iVector<vtype,Ns> &ret, const iVector<vtype,Ns> &rhs){
ret(0) =-rhs(0);
ret(1) =-rhs(1);
ret(2) = rhs(2);
ret(3) = rhs(3);
};
template<class vtype, int N>
inline auto operator*(const GammaL &gl, const iMatrix<vtype, N> &arg)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=gl*arg._internal[i][j];
}}
return ret;
}
//general right multiply
template<class vtype>
inline auto operator*(const iScalar<vtype> &arg, const GammaL &gl)
->typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype>>::type
{
iScalar<vtype> ret;
ret._internal=arg._internal*gl;
return ret;
}
template<class vtype, int N>
inline auto operator * (const iMatrix<vtype, N> &arg, const GammaL &gl)
->typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N>>::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=arg._internal[i][j]*gl;
}}
return ret;
}
#ifdef GRID_WARN_SUBOPTIMAL
#warning "Optimisation alert switch over to multGammaX early "
#endif
}}
///////////////////////////////////////////////////////////////////////////////////////////////////
// Operator * : first case this is not a spin index, so recurse
///////////////////////////////////////////////////////////////////////////////////////////////////
// FIXME
//
// Optimisation; switch over to a "multGammaX(ret._internal,arg._internal)" style early and
// note that doing so from the lattice operator will avoid copy back and case switch overhead, as
// was done for the tensor math operator to remove operator * notation early
//
//left multiply
template<class vtype> inline auto operator * ( const Gamma &G,const iScalar<vtype> &arg) ->
typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type
{
iScalar<vtype> ret;
ret._internal=G*arg._internal;
return ret;
}
template<class vtype,int N> inline auto operator * ( const Gamma &G,const iVector<vtype,N> &arg) ->
typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal[i]=G*arg._internal[i];
}
return ret;
}
template<class vtype,int N> inline auto operator * ( const Gamma &G,const iMatrix<vtype,N> &arg) ->
typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=G*arg._internal[i][j];
}}
return ret;
}
//right multiply
template<class vtype> inline auto operator * (const iScalar<vtype> &arg, const Gamma &G) ->
typename std::enable_if<matchGridTensorIndex<iScalar<vtype>,SpinorIndex>::notvalue,iScalar<vtype> >::type
{
iScalar<vtype> ret;
ret._internal=arg._internal*G;
return ret;
}
template<class vtype,int N> inline auto operator * (const iVector<vtype,N> &arg, const Gamma &G) ->
typename std::enable_if<matchGridTensorIndex<iVector<vtype,N>,SpinorIndex>::notvalue,iVector<vtype,N> >::type
{
iVector<vtype,N> ret;
for(int i=0;i<N;i++){
ret._internal=arg._internal[i]*G;
}
return ret;
}
template<class vtype,int N> inline auto operator * (const iMatrix<vtype,N> &arg, const Gamma &G) ->
typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,N>,SpinorIndex>::notvalue,iMatrix<vtype,N> >::type
{
iMatrix<vtype,N> ret;
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
ret._internal[i][j]=arg._internal[i][j]*G;
}}
return ret;
}
////////////////////////////////////////////////////////
// When we hit the spin index this matches and we stop
////////////////////////////////////////////////////////
template<class vtype> inline auto operator * ( const Gamma &G,const iMatrix<vtype,Ns> &arg) ->
typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type
{
iMatrix<vtype,Ns> ret;
switch (G._g) {
case Gamma::Identity:
ret = arg;
break;
case Gamma::MinusIdentity:
ret = -arg;
break;
case Gamma::GammaX:
multGammaX(ret,arg);
break;
case Gamma::MinusGammaX:
multMinusGammaX(ret,arg);
break;
case Gamma::GammaY:
multGammaY(ret,arg);
break;
case Gamma::MinusGammaY:
multMinusGammaY(ret,arg);
break;
case Gamma::GammaZ:
multGammaZ(ret,arg);
break;
case Gamma::MinusGammaZ:
multMinusGammaZ(ret,arg);
break;
case Gamma::GammaT:
multGammaT(ret,arg);
break;
case Gamma::MinusGammaT:
multMinusGammaT(ret,arg);
break;
case Gamma::Gamma5:
multGamma5(ret,arg);
break;
case Gamma::MinusGamma5:
multMinusGamma5(ret,arg);
break;
default:
assert(0);
break;
}
return ret;
}
// Could have used type trait for Matrix/vector and then an enable if to share code
template<class vtype> inline auto operator * ( const Gamma &G,const iVector<vtype,Ns> &arg) ->
typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type
{
iVector<vtype,Ns> ret;
switch (G._g) {
case Gamma::Identity:
ret = arg;
break;
case Gamma::MinusIdentity:
ret = -arg;
break;
case Gamma::GammaX:
multGammaX(ret,arg);
break;
case Gamma::MinusGammaX:
multMinusGammaX(ret,arg);
break;
case Gamma::GammaY:
multGammaY(ret,arg);
break;
case Gamma::MinusGammaY:
multMinusGammaY(ret,arg);
break;
case Gamma::GammaZ:
multGammaZ(ret,arg);
break;
case Gamma::MinusGammaZ:
multMinusGammaZ(ret,arg);
break;
case Gamma::GammaT:
multGammaT(ret,arg);
break;
case Gamma::MinusGammaT:
multMinusGammaT(ret,arg);
break;
case Gamma::Gamma5:
multGamma5(ret,arg);
break;
case Gamma::MinusGamma5:
multMinusGamma5(ret,arg);
break;
default:
assert(0);
break;
}
return ret;
}
template<class vtype> inline auto operator * (const iMatrix<vtype,Ns> &arg, const Gamma &G) ->
typename std::enable_if<matchGridTensorIndex<iMatrix<vtype,Ns>,SpinorIndex>::value,iMatrix<vtype,Ns> >::type
{
iMatrix<vtype,Ns> ret;
switch (G._g) {
case Gamma::Identity:
ret = arg;
break;
case Gamma::MinusIdentity:
ret = -arg;
break;
case Gamma::GammaX:
rmultGammaX(ret,arg);
break;
case Gamma::MinusGammaX:
rmultMinusGammaX(ret,arg);
break;
case Gamma::GammaY:
rmultGammaY(ret,arg);
break;
case Gamma::MinusGammaY:
rmultMinusGammaY(ret,arg);
break;
case Gamma::GammaZ:
rmultGammaZ(ret,arg);
break;
case Gamma::MinusGammaZ:
rmultMinusGammaZ(ret,arg);
break;
case Gamma::GammaT:
rmultGammaT(ret,arg);
break;
case Gamma::MinusGammaT:
rmultMinusGammaT(ret,arg);
break;
case Gamma::Gamma5:
rmultGamma5(ret,arg);
break;
case Gamma::MinusGamma5:
rmultMinusGamma5(ret,arg);
break;
default:
assert(0);
break;
}
return ret;
}
/* Output from test
./Grid_gamma
Identity((1,0),(0,0),(0,0),(0,0))
((0,0),(1,0),(0,0),(0,0))
((0,0),(0,0),(1,0),(0,0))
((0,0),(0,0),(0,0),(1,0)) OK
GammaX ((0,0),(0,0),(0,0),(0,1))
((0,0),(0,0),(0,1),(0,0))
((0,0),(0,-1),(0,0),(0,0))
((0,-1),(0,0),(0,0),(0,0)) OK
* Gx
* 0 0 0 i
* 0 0 i 0
* 0 -i 0 0
* -i 0 0 0
GammaY ((-0,-0),(-0,-0),(-0,-0),(-1,-0))
((0,0),(0,0),(1,0),(0,0))
((0,0),(1,0),(0,0),(0,0)) OK
((-1,-0),(-0,-0),(-0,-0),(-0,-0))
*Gy
* 0 0 0 -1 [0] -+ [3]
* 0 0 1 0 [1] +- [2]
* 0 1 0 0
* -1 0 0 0
GammaZ ((0,0),(0,0),(0,1),(0,0))
((0,0),(0,0),(0,0),(0,-1))
((0,-1),(0,0),(0,0),(0,0))
((0,0),(0,1),(0,0),(0,0)) OK
* 0 0 i 0 [0]+-i[2]
* 0 0 0 -i [1]-+i[3]
* -i 0 0 0
* 0 i 0 0
GammaT ((0,0),(0,0),(1,0),(0,0))
((0,0),(0,0),(0,0),(1,0)) OK
((1,0),(0,0),(0,0),(0,0))
((0,0),(1,0),(0,0),(0,0))
* 0 0 1 0 [0]+-[2]
* 0 0 0 1 [1]+-[3]
* 1 0 0 0
* 0 1 0 0
Gamma5 ((1,0),(0,0),(0,0),(0,0))
((0,0),(1,0),(0,0),(0,0))
((-0,-0),(-0,-0),(-1,-0),(-0,-0))
((-0,-0),(-0,-0),(-0,-0),(-1,-0))
* 1 0 0 0 [0]+-[2]
* 0 1 0 0 [1]+-[3] OK
* 0 0 -1 0
* 0 0 0 -1
*/
} //namespace QCD
} // Grid
#endif

1140
lib/qcd/spin/Gamma.cc Normal file
View File

File diff suppressed because it is too large Load Diff

1348
lib/qcd/spin/Gamma.h Normal file
View File

File diff suppressed because it is too large Load Diff

1488
lib/qcd/spin/gamma-gen/gamma-gen.nb Normal file
View File

File diff suppressed because it is too large Load Diff

10
lib/qcd/utils/LinalgUtils.h
View File

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

2
lib/qcd/utils/SpaceTimeGrid.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
namespace Grid {
namespace QCD {

2
lib/serialisation/BinaryIO.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
using namespace Grid;
using namespace std;

2
lib/serialisation/Hdf5IO.cc
View File

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

17
lib/serialisation/MacroMagic.h
View File

@@ -150,14 +150,14 @@ friend inline bool operator==(const cname &lhs, const cname &rhs) {\
class name: public Grid::Serializable\
{\
public:\
enum EnumType\
enum\
{\
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMVAL,__VA_ARGS__))\
undefname = -1\
};\
public:\
name(void): value_(undefname) {};\
name(EnumType value): value_(value) {};\
name(int value): value_(value) {};\
template <typename T>\
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__))\
else {obj = name::undefname;}\
}\
inline operator EnumType(void) const\
inline operator int(void) const\
{\
return value_;\
}\
@@ -189,8 +189,17 @@ public:\
}\
return os;\
}\
inline friend std::istream & operator>>(std::istream &is, name &obj)\
{\
std::string buf;\
is >> buf;\
if (buf == #undefname) {obj = name::undefname;}\
GRID_MACRO_EVAL(GRID_MACRO_MAP(GRID_MACRO_ENUMTEST,__VA_ARGS__))\
else {obj = name::undefname;}\
return is;\
}\
private:\
EnumType value_;\
int value_;\
};

3
lib/serialisation/TextIO.cc
View File

@@ -27,8 +27,7 @@
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <ctype.h>
#include <Grid/Grid.h>
using namespace Grid;
using namespace std;

2
lib/serialisation/XmlIO.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
using namespace Grid;
using namespace std;

390
lib/sitmo_rng/sitmo_prng_engine.hpp Normal file
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 e’s 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 e’s 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

2
lib/stencil/Lebesgue.cc
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
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <Grid/Grid.h>
#include <algorithm>
namespace Grid {

2
lib/stencil/Stencil_common.cc
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
*************************************************************************************/
/* END LEGAL */
#include "Grid.h"
#include <Grid/Grid.h>
namespace Grid {
}