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

Added generators for the adjoint representation

Browse Source
This commit is contained in:
Guido Cossu 2016-07-08 15:40:11 +01:00
parent fbf96b1bbb
commit 5028969d4b
15 changed files with 780 additions and 665 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
benchmarks/Benchmark_dwf_ntpf
View File

Binary file not shown.

BIN
benchmarks/Benchmark_zmm
View File

Binary file not shown.

72
lib/qcd/hmc/integrators/Integrator.h
View File

@ -25,7 +25,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
//--------------------------------------------------------------------
@ -48,31 +49,24 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
namespace QCD {
struct IntegratorParameters {
int Nexp;
int MDsteps; // number of outer steps
RealD trajL; // trajectory length
RealD stepsize;
IntegratorParameters(int MDsteps_,
RealD trajL_=1.0,
int Nexp_=12):
Nexp(Nexp_),
IntegratorParameters(int MDsteps_, RealD trajL_ = 1.0, int Nexp_ = 12)
: Nexp(Nexp_),
MDsteps(MDsteps_),
trajL(trajL_),
stepsize(trajL/MDsteps)
{
stepsize(trajL / MDsteps){
// empty body constructor
};
};
/*! @brief Class for Molecular Dynamics management */
template <class GaugeField, class SmearingPolicy>
class Integrator {
protected:
typedef IntegratorParameters ParameterType;
IntegratorParameters Params;
@ -89,17 +83,18 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
// Should match any legal (SU(n)) gauge field
// Need to use this template to match Ncol to pass to SU<N> class
template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){
template <int Ncol, class vec>
void generate_momenta(Lattice<iVector<iScalar<iMatrix<vec, Ncol> >, Nd> >& P,
GridParallelRNG& pRNG) {
typedef Lattice<iScalar<iScalar<iMatrix<vec, Ncol> > > > GaugeLinkField;
GaugeLinkField Pmu(P._grid);
Pmu = zero;
for (int mu = 0; mu < Nd; mu++) {
SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
SU<Ncol>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
PokeIndex<LorentzIndex>(P, Pmu, mu);
}
}
// ObserverList observers; // not yet
// typedef std::vector<Observer*> ObserverList;
// void register_observers();
@ -109,7 +104,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
t_P[level] += ep;
update_P(P, U, level, ep);
std::cout<<GridLogIntegrator<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
std::cout << GridLogIntegrator << "[" << level << "] P "
<< " dt " << ep << " : t_P " << t_P[level] << std::endl;
}
void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) {
@ -119,10 +115,14 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta
std::cout<< GridLogIntegrator << "Smearing (on/off): "<<as[level].actions.at(a)->is_smeared <<std::endl;
std::cout << GridLogIntegrator
<< "Smearing (on/off): " << as[level].actions.at(a)->is_smeared
<< std::endl;
if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
force = Ta(force);
std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <<std::endl;
std::cout << GridLogIntegrator
<< "Force average: " << norm2(force) / (U._grid->gSites())
<< std::endl;
Mom -= force * ep;
}
}
@ -132,8 +132,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
t_U += ep;
int fl = levels - 1;
std::cout<< GridLogIntegrator <<" "<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
std::cout << GridLogIntegrator << " "
<< "[" << fl << "] U "
<< " dt " << ep << " : t_U " << t_U << std::endl;
}
void update_U(GaugeField& Mom, GaugeField& U, double ep) {
// rewrite exponential to deal automatically with the lorentz index?
@ -153,17 +154,9 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
virtual void step(GaugeField& U, int level, int first, int last) = 0;
public:
Integrator(GridBase* grid,
IntegratorParameters Par,
ActionSet<GaugeField> & Aset,
SmearingPolicy &Sm):
Params(Par),
as(Aset),
P(grid),
levels(Aset.size()),
Smearer(Sm)
{
Integrator(GridBase* grid, IntegratorParameters Par,
ActionSet<GaugeField>& Aset, SmearingPolicy& Sm)
: Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm) {
t_P.resize(levels, 0.0);
t_U = 0.0;
// initialization of smearer delegated outside of Integrator
@ -179,7 +172,8 @@ public:
for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
// get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID
GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
GaugeField& Us =
Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
as[level].actions.at(actionID)->refresh(Us, pRNG);
}
}
@ -188,7 +182,8 @@ public:
// Calculate action
RealD S(GaugeField& U) { // here also U not used
LatticeComplex Hloc(U._grid); Hloc = zero;
LatticeComplex Hloc(U._grid);
Hloc = zero;
// Momenta
for (int mu = 0; mu < Nd; mu++) {
auto Pmu = PeekIndex<LorentzIndex>(P, mu);
@ -205,9 +200,11 @@ public:
for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
// get gauge field from the SmearingPolicy and
// based on the boolean is_smeared in actionID
GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
GaugeField& Us =
Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
Hterm = as[level].actions.at(actionID)->S(Us);
std::cout<<GridLogMessage << "S Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
std::cout << GridLogMessage << "S Level " << level << " term "
<< actionID << " H = " << Hterm << std::endl;
H += Hterm;
}
}
@ -216,7 +213,6 @@ public:
}
void integrate(GaugeField& U) {
// reset the clocks
t_U = 0;
for (int level = 0; level < as.size(); ++level) {
@ -232,16 +228,14 @@ public:
// Check the clocks all match on all levels
for (int level = 0; level < as.size(); ++level) {
assert(fabs(t_U - t_P[level]) < 1.0e-6); // must be the same
std::cout<<GridLogIntegrator<<" times["<<level<<"]= "<<t_P[level]<< " " << t_U <<std::endl;
std::cout << GridLogIntegrator << " times[" << level
<< "]= " << t_P[level] << " " << t_U << std::endl;
}
// and that we indeed got to the end of the trajectory
assert(fabs(t_U - Params.trajL) < 1.0e-6);
}
};
}
}
#endif // INTEGRATOR_INCLUDED

298
lib/qcd/utils/SUn.h
View File

@ -25,7 +25,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
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
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#ifndef QCD_UTIL_SUN_H
@ -37,15 +38,19 @@ namespace Grid {
template <int ncolour>
class SU {
public:
static int generators(void) { return ncolour * ncolour - 1; }
static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; }
template<typename vtype> using iSUnMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > > ;
template<typename vtype> using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > > ;
template <typename vtype>
using iSUnMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
template <typename vtype>
using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
template <typename vtype>
using iSUnAdjointMatrix = iScalar<iScalar<iMatrix<vtype, (ncolour*ncolour - 1)> > >;
//////////////////////////////////////////////////////////////////////////////////////////////////
// Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, SU<2>::LatticeMatrix etc...
// Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
// SU<2>::LatticeMatrix etc...
//////////////////////////////////////////////////////////////////////////////////////////////////
typedef iSUnMatrix<Complex> Matrix;
typedef iSUnMatrix<ComplexF> MatrixF;
@ -55,6 +60,16 @@ public:
typedef iSUnMatrix<vComplexF> vMatrixF;
typedef iSUnMatrix<vComplexD> vMatrixD;
// Actually the adjoint matrices are real...
// Consider this overhead... FIXME
typedef iSUnAdjointMatrix<Complex> AMatrix;
typedef iSUnAdjointMatrix<ComplexF> AMatrixF;
typedef iSUnAdjointMatrix<ComplexD> AMatrixD;
typedef iSUnAdjointMatrix<vComplex> vAMatrix;
typedef iSUnAdjointMatrix<vComplexF> vAMatrixF;
typedef iSUnAdjointMatrix<vComplexD> vAMatrixD;
typedef Lattice<vMatrix> LatticeMatrix;
typedef Lattice<vMatrixF> LatticeMatrixF;
typedef Lattice<vMatrixD> LatticeMatrixD;
@ -71,13 +86,13 @@ public:
typedef Lattice<vSU2MatrixF> LatticeSU2MatrixF;
typedef Lattice<vSU2MatrixD> LatticeSU2MatrixD;
////////////////////////////////////////////////////////////////////////
// There are N^2-1 generators for SU(N).
//
// We take a traceless hermitian generator basis as follows
//
// * Normalisation: trace ta tb = 1/2 delta_ab
// * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab
// T_F = 1/2 for SU(N) groups
//
// * Off diagonal
// - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y
@ -93,11 +108,14 @@ public:
// - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc
//
// - We enumerate the row-col pairs.
// - for each row col pair there is a (sigma_x) and a (sigma_y) like generator
// - for each row col pair there is a (sigma_x) and a (sigma_y) like
// generator
//
//
// t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => delta_{i,i1} delta_{j,i2} + delta_{i,i1} delta_{j,i2}
// t^a_ij = { in Nc(Nc-1)/2 ... Nc^(Nc-1) -1} => i delta_{i,i1} delta_{j,i2} - i delta_{i,i1} delta_{j,i2}
// t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} +
// delta_{i,i1} delta_{j,i2})
// t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1}
// delta_{j,i2} - i delta_{i,i1} delta_{j,i2})
//
// * Diagonal; must be traceless and normalised
// - Sequence is
@ -114,6 +132,20 @@ public:
// ( 1 )
// ( 1 ) / sqrt(3) /2 = 1/2 lambda_8
// ( -2)
//
//
// * Adjoint representation generators
//
// base for NxN hermitian traceless matrices
// normalized to 1:
//
// (e_Adj)^a = t^a / sqrt(T_F)
//
// then the real, antisymmetric generators for the adjoint representations
// are computed ( shortcut: e^a == (e_Adj)^a )
//
// (iT_adj)^d_ab = i tr[e^a t^d e^b - t^d e^a e^b]
//
////////////////////////////////////////////////////////////////////////
template <class cplx>
static void generator(int lieIndex, iSUnMatrix<cplx> &ta) {
@ -127,10 +159,12 @@ public:
generatorDiagonal(diagIndex, ta);
return;
}
sigxy = lieIndex&0x1;
sigxy = lieIndex & 0x1;//even or odd
su2Index = lieIndex >> 1;
if ( sigxy ) generatorSigmaY(su2Index,ta);
else generatorSigmaX(su2Index,ta);
if (sigxy)
generatorSigmaY(su2Index, ta);
else
generatorSigmaX(su2Index, ta);
}
template <class cplx>
static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
@ -153,24 +187,45 @@ public:
}
template <class cplx>
static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) {
// diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...)
ta = zero;
int trsq=0;
int last=diagIndex+1;
for(int i=0;i<=diagIndex;i++){
int k = diagIndex + 1;// diagIndex starts from 0
for (int i = 0; i <= diagIndex; i++) {// k iterations
ta()()(i, i) = 1.0;
trsq++;
}
ta()()(last,last) = -last;
trsq+=last*last;
RealD nrm = 1.0/std::sqrt(2.0*trsq);
ta()()(k,k) = -k;//indexing starts from 0
RealD nrm = 1.0 / std::sqrt(2.0 * k*(k+1));
ta = ta * nrm;
}
template <class cplx>
static void generatorAdjoint(int Index, iSUnAdjointMatrix<cplx> &iAdjTa){
// returns i(T_Adj)^index necessary for the projectors
// see definitions above
iAdjTa = zero;
Vector< iSUnMatrix<cplx> > ta(ncolour*ncolour -1);
iSUnMatrix<cplx> tmp;
// FIXME not very efficient to get all the generators everytime
for (int a = 0; a < (ncolour * ncolour - 1); a++)
generator(a, ta[a]);
for (int a = 0; a < (ncolour*ncolour - 1); a++){
tmp = ta[a] * ta[Index] - ta[Index] * ta[a];
for (int b = 0; b < (ncolour*ncolour - 1); b++){
iSUnMatrix<cplx> tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization
Complex iTr = TensorRemove(timesI(trace(tmp1)));
iAdjTa()()(b,a) = iTr;
}
}
}
////////////////////////////////////////////////////////////////////////
// Map a su2 subgroup number to the pair of rows that are non zero
////////////////////////////////////////////////////////////////////////
static void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2));
int spare = su2_index;
@ -187,8 +242,7 @@ public:
static void su2Extract(Lattice<iSinglet<vcplx> > &Determinant,
Lattice<iSU2Matrix<vcplx> > &subgroup,
const Lattice<iSUnMatrix<vcplx> > &source,
int su2_index)
{
int su2_index) {
GridBase *grid(source._grid);
conformable(subgroup, source);
conformable(subgroup, Determinant);
@ -209,11 +263,9 @@ PARALLEL_FOR_LOOP
subgroup._odata[ss] = Sigma;
// this should be purely real
Determinant._odata[ss] = Sigma()()(0,0)*Sigma()()(1,1)
- Sigma()()(0,1)*Sigma()()(1,0);
Determinant._odata[ss] =
Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0);
}
}
//////////////////////////////////////////////////////////////////////////////////////////
@ -221,9 +273,7 @@ PARALLEL_FOR_LOOP
//////////////////////////////////////////////////////////////////////////////////////////
template <class vcplx>
static void su2Insert(const Lattice<iSU2Matrix<vcplx> > &subgroup,
Lattice<iSUnMatrix<vcplx> > &dest,
int su2_index)
{
Lattice<iSUnMatrix<vcplx> > &dest, int su2_index) {
GridBase *grid(dest._grid);
conformable(subgroup, dest);
int i0, i1;
@ -239,24 +289,22 @@ PARALLEL_FOR_LOOP
}
}
///////////////////////////////////////////////
// Generate e^{ Re Tr Staple Link} dlink
//
// *** Note Staple should be appropriate linear compbination between all staples.
// *** Note Staple should be appropriate linear compbination between all
// staples.
// *** If already by beta pass coefficient 1.0.
// *** This routine applies the additional 1/Nc factor that comes after trace in action.
// *** This routine applies the additional 1/Nc factor that comes after trace
// in action.
//
///////////////////////////////////////////////
static void SubGroupHeatBath( GridSerialRNG &sRNG,
GridParallelRNG &pRNG,
static void SubGroupHeatBath(
GridSerialRNG &sRNG, GridParallelRNG &pRNG,
RealD beta, // coeff multiplying staple in action (with no 1/Nc)
LatticeMatrix &link,
const LatticeMatrix &barestaple, // multiplied by action coeffs so th
int su2_subgroup,
int nheatbath,
LatticeInteger &wheremask)
{
int su2_subgroup, int nheatbath, LatticeInteger &wheremask) {
GridBase *grid = link._grid;
int ntrials = 0;
@ -267,22 +315,30 @@ PARALLEL_FOR_LOOP
staple = barestaple * (beta / ncolour);
LatticeMatrix V(grid); V = link*staple;
LatticeMatrix V(grid);
V = link * staple;
// Subgroup manipulation in the lie algebra space
LatticeSU2Matrix u(grid); // Kennedy pendleton "u" real projected normalised Sigma
LatticeSU2Matrix u(
grid); // Kennedy pendleton "u" real projected normalised Sigma
LatticeSU2Matrix uinv(grid);
LatticeSU2Matrix ua(grid); // a in pauli form
LatticeSU2Matrix b(grid); // rotated matrix after hb
// Some handy constant fields
LatticeComplex ones (grid); ones = 1.0;
LatticeComplex zeros(grid); zeros=zero;
LatticeReal rones (grid); rones = 1.0;
LatticeReal rzeros(grid); rzeros=zero;
LatticeComplex ones(grid);
ones = 1.0;
LatticeComplex zeros(grid);
zeros = zero;
LatticeReal rones(grid);
rones = 1.0;
LatticeReal rzeros(grid);
rzeros = zero;
LatticeComplex udet(grid); // determinant of real(staple)
LatticeInteger mask_true (grid); mask_true = 1;
LatticeInteger mask_false(grid); mask_false= 0;
LatticeInteger mask_true(grid);
mask_true = 1;
LatticeInteger mask_false(grid);
mask_false = 0;
/*
PLB 156 P393 (1985) (Kennedy and Pendleton)
@ -299,7 +355,8 @@ Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' "
beta S = const - beta/Nc Re Tr h Sigma'
= const - Re Tr h Sigma
Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex arbitrary.
Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex
arbitrary.
Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij
Re Tr h Sigma = 2 h_j Re Sigma_j
@ -335,9 +392,12 @@ Note: Product b' xi is unvariant because scaling Sigma leaves
LatticeSU2Matrix lident(grid);
SU2Matrix ident = Complex(1.0);
SU2Matrix pauli1; SU<2>::generator(0,pauli1);
SU2Matrix pauli2; SU<2>::generator(1,pauli2);
SU2Matrix pauli3; SU<2>::generator(2,pauli3);
SU2Matrix pauli1;
SU<2>::generator(0, pauli1);
SU2Matrix pauli2;
SU<2>::generator(1, pauli2);
SU2Matrix pauli3;
SU<2>::generator(2, pauli3);
pauli1 = timesI(pauli1) * 2.0;
pauli2 = timesI(pauli2) * 2.0;
pauli3 = timesI(pauli3) * 2.0;
@ -352,7 +412,8 @@ Note: Product b' xi is unvariant because scaling Sigma leaves
udet = where(adet > machine_epsilon, udet, cone);
xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag]
u = 0.5*u*pow(xi,-1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag]
u = 0.5 * u *
pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag]
// Debug test for sanity
uinv = adj(u);
@ -363,19 +424,24 @@ Note: Product b' xi is unvariant because scaling Sigma leaves
Measure: Haar measure dh has d^4a delta(1-|a^2|)
In polars:
da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2)
= da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + r) )
= da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) +
r) )
= da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) )
Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters through xi
Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters
through xi
= e^{2 xi (h.u)} dh
= e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi h2u2}.e^{2 xi h3u3} dh
= e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi
h2u2}.e^{2 xi h3u3} dh
Therefore for each site, take xi for that site
i) generate |a0|<1 with dist
(1-a0^2)^0.5 e^{2 xi a0 } da0
Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc factor in Chroma ]
A. Generate two uniformly distributed pseudo-random numbers R and R', R'', R''' in the unit interval;
Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc
factor in Chroma ]
A. Generate two uniformly distributed pseudo-random numbers R and R', R'',
R''' in the unit interval;
B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha;
C. Set C = cos^2(2pi R"), with R" another uniform random number in [0,1] ;
D. Set A = XC;
@ -383,7 +449,8 @@ E. Let d = X'+A;
F. If R'''^2 :> 1 - 0.5 d, go back to A;
G. Set a0 = 1 - d;
Note that in step D setting B ~ X - A and using B in place of A in step E will generate a second independent a 0 value.
Note that in step D setting B ~ X - A and using B in place of A in step E will
generate a second independent a 0 value.
*/
/////////////////////////////////////////////////////////
@ -394,20 +461,22 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
rtmp = where(wheremask, rones, rzeros);
RealD numSites = sum(rtmp);
RealD numAccepted;
LatticeInteger Accepted(grid); Accepted = zero;
LatticeInteger Accepted(grid);
Accepted = zero;
LatticeInteger newlyAccepted(grid);
std::vector<LatticeReal> xr(4, grid);
std::vector<LatticeReal> a(4, grid);
LatticeReal d(grid); d=zero;
LatticeReal d(grid);
d = zero;
LatticeReal alpha(grid);
// std::cout<<GridLogMessage<<"xi "<<xi <<std::endl;
alpha = toReal(2.0 * xi);
do {
// A. Generate two uniformly distributed pseudo-random numbers R and R', R'', R''' in the unit interval;
// A. Generate two uniformly distributed pseudo-random numbers R and R',
// R'', R''' in the unit interval;
random(pRNG, xr[0]);
random(pRNG, xr[1]);
random(pRNG, xr[2]);
@ -430,10 +499,13 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
d = where(Accepted, d, xr[2] + xr[1] * xr[3]);
// F. If R'''^2 :> 1 - 0.5 d, go back to A;
LatticeReal thresh(grid); thresh = 1.0-d*0.5;
LatticeReal thresh(grid);
thresh = 1.0 - d * 0.5;
xrsq = xr[0] * xr[0];
LatticeInteger ione(grid); ione = 1;
LatticeInteger izero(grid); izero=zero;
LatticeInteger ione(grid);
ione = 1;
LatticeInteger izero(grid);
izero = zero;
newlyAccepted = where(xrsq < thresh, ione, izero);
Accepted = where(newlyAccepted, newlyAccepted, Accepted);
@ -462,18 +534,18 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeReal sin_theta(grid);
LatticeReal phi(grid);
random(pRNG,phi); phi = phi * twopi; // uniform in [0,2pi]
random(pRNG,cos_theta); cos_theta=(cos_theta*2.0)-1.0; // uniform in [-1,1]
random(pRNG, phi);
phi = phi * twopi; // uniform in [0,2pi]
random(pRNG, cos_theta);
cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1]
sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta));
a[1] = a123mag * sin_theta * cos(phi);
a[2] = a123mag * sin_theta * sin(phi);
a[3] = a123mag * cos_theta;
ua = toComplex(a[0])*ident
+ toComplex(a[1])*pauli1
+ toComplex(a[2])*pauli2
+ toComplex(a[3])*pauli3;
ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 +
toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3;
b = 1.0;
b = where(wheremask, uinv * ua, b);
@ -508,52 +580,87 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
/////////////////////////////////
}
static void printGenerators(void)
{
static void printGenerators(void) {
for (int gen = 0; gen < generators(); gen++) {
Matrix ta;
generator(gen, ta);
std::cout<<GridLogMessage<< "Nc = "<<ncolour<<" t_"<<gen<<std::endl;
std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
<< std::endl;
std::cout << GridLogMessage << ta << std::endl;
}
}
static void printAdjointGenerators(void) {
for (int gen = 0; gen < generators(); gen++) {
AMatrix ta;
generatorAdjoint(gen, ta);
std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
<< std::endl;
std::cout << GridLogMessage << ta << std::endl;
}
}
static void testGenerators(void) {
Matrix ta;
Matrix tb;
std::cout<<GridLogMessage<<"Checking trace ta tb is 0.5 delta_ab"<<std::endl;
std::cout << GridLogMessage << "Fundamental - Checking trace ta tb is 0.5 delta_ab"
<< std::endl;
for (int a = 0; a < generators(); a++) {
for (int b = 0; b < generators(); b++) {
generator(a, ta);
generator(b, tb);
Complex tr = TensorRemove(trace(ta * tb));
std::cout<<GridLogMessage<<tr<<" ";
std::cout << GridLogMessage << "("<< a << "," << b << ") = "<< tr << std::endl;
if (a == b) assert(abs(tr - Complex(0.5)) < 1.0e-6);
if (a != b) assert(abs(tr) < 1.0e-6);
}
std::cout << GridLogMessage << std::endl;
}
std::cout<<GridLogMessage<<"Checking hermitian"<<std::endl;
std::cout << GridLogMessage << "Fundamental - Checking if hermitian" << std::endl;
for (int a = 0; a < generators(); a++) {
generator(a, ta);
std::cout<<GridLogMessage<<a<<" ";
std::cout << GridLogMessage << a << std::endl;
assert(norm2(ta - adj(ta)) < 1.0e-6) ;
}
std::cout << GridLogMessage << std::endl;
std::cout<<GridLogMessage<<"Checking traceless"<<std::endl;
std::cout << GridLogMessage << "Fundamental - Checking if traceless" << std::endl;
for (int a = 0; a < generators(); a++) {
generator(a, ta);
Complex tr = TensorRemove(trace(ta));
std::cout<<GridLogMessage<<a<<" ";
std::cout << GridLogMessage << a << " " << std::endl;
assert(abs(tr) < 1.0e-6);
}
std::cout << GridLogMessage << std::endl;
AMatrix adjTa;
std::cout << GridLogMessage << "Adjoint - Checking if real" << std::endl;
for (int a = 0; a < generators(); a++) {
generatorAdjoint(a, adjTa);
std::cout << GridLogMessage << a << std::endl;
assert(norm2(adjTa - conjugate(adjTa)) < 1.0e-6);
}
std::cout << GridLogMessage << std::endl;
std::cout << GridLogMessage << "Adjoint - Checking if antisymmetric" << std::endl;
for (int a = 0; a < generators(); a++) {
generatorAdjoint(a, adjTa);
std::cout << GridLogMessage << a << std::endl;
assert(norm2(adjTa + transpose(adjTa)) < 1.0e-6);
}
std::cout << GridLogMessage << std::endl;
}
// reunitarise??
template <typename LatticeMatrixType>
static void LieRandomize(GridParallelRNG &pRNG,LatticeMatrixType &out,double scale=1.0){
static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out,
double scale = 1.0) {
GridBase *grid = out._grid;
typedef typename LatticeMatrixType::vector_type vector_type;
@ -562,7 +669,8 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
typedef iSinglet<vector_type> vTComplexType;
typedef Lattice<vTComplexType> LatticeComplexType;
typedef typename GridTypeMapper<typename LatticeMatrixType::vector_object>::scalar_object MatrixType;
typedef typename GridTypeMapper<
typename LatticeMatrixType::vector_object>::scalar_object MatrixType;
LatticeComplexType ca(grid);
LatticeMatrixType lie(grid);
@ -573,7 +681,6 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
lie = zero;
for (int a = 0; a < generators(); a++) {
random(pRNG, ca);
ca = (ca + conjugate(ca)) * 0.5;
@ -588,7 +695,9 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
taExp(lie, out);
}
static void GaussianLieAlgebraMatrix(GridParallelRNG &pRNG,LatticeMatrix &out,double scale=1.0){
static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG,
LatticeMatrix &out,
Real scale = 1.0) {
GridBase *grid = out._grid;
LatticeReal ca(grid);
LatticeMatrix la(grid);
@ -602,7 +711,20 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
la = toComplex(ca) * ci * ta;
out += la;
}
}
static void FundamentalLieAlgebraMatrix(Vector<Real> &h, LatticeMatrix &out,
Real scale = 1.0) {
GridBase *grid = out._grid;
LatticeMatrix la(grid);
Matrix ta;
out = zero;
for (int a = 0; a < generators(); a++) {
generator(a, ta);
la = Complex(0.0, h[a]) * scale * ta;
out += la;
}
}
template <typename GaugeField>
@ -617,7 +739,8 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
PokeIndex<LorentzIndex>(out, Umu, mu);
}
}
static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){
static void TepidConfiguration(GridParallelRNG &pRNG,
LatticeGaugeField &out) {
LatticeMatrix Umu(out._grid);
for (int mu = 0; mu < Nd; mu++) {
LieRandomize(pRNG, Umu, 0.01);
@ -646,21 +769,18 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
ex = xn + ComplexType(1.0); // 1+x
// Do a 12th order exponentiation
for(int i=2; i <= 12; ++i)
{
for (int i = 2; i <= 12; ++i) {
nfac = nfac / RealD(i); // 1/2, 1/2.3 ...
xn = xn * x; // x2, x3,x4....
ex = ex + xn * nfac; // x2/2!, x3/3!....
}
}
};
typedef SU<2> SU2;
typedef SU<3> SU3;
typedef SU<4> SU4;
typedef SU<5> SU5;
}
}
#endif

2
tests/Test_contfrac_force.cc
View File

@ -96,7 +96,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

2
tests/Test_dwf_force.cc
View File

@ -96,7 +96,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

2
tests/Test_dwf_gpforce.cc
View File

@ -113,7 +113,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
Hmom -= real(sum(trace(mommu*mommu)));

2
tests/Test_gp_rect_force.cc
View File

@ -82,7 +82,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/Test_gpdwf_force.cc
View File

@ -100,7 +100,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
@ -169,7 +169,7 @@ int main (int argc, char ** argv)
//
// Pmu = zero;
// for(int mu=0;mu<Nd;mu++){
// SU<Ncol>::GaussianLieAlgebraMatrix(pRNG, Pmu);
// SU<Ncol>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
// PokeIndex<LorentzIndex>(P, Pmu, mu);
// }
//

41
tests/Test_lie_generators.cc
View File

@ -23,48 +23,52 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
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
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid.h>
#include <qcd/utils/CovariantCshift.h>
#include <qcd/utils/WilsonLoops.h>
#include <qcd/utils/SUn.h>
#include <qcd/utils/WilsonLoops.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
int main(int argc, char** argv) {
Grid_init(&argc, &argv);
std::vector<int> latt({4, 4, 4, 8});
GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplex::Nsimd()),
GridDefaultMpi());
GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for SU(2)" << std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
SU2::printGenerators();
SU2::printAdjointGenerators();
SU2::testGenerators();
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for SU(3)" << std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
SU3::printGenerators();
SU3::printAdjointGenerators();
SU3::testGenerators();
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"* Generators for SU(4)"<<std::endl;
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// SU4::printGenerators();
// SU4::testGenerators();
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Generators for SU(4)"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
SU4::printGenerators();
SU4::testGenerators();
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"* Generators for SU(5)"<<std::endl;
@ -72,8 +76,5 @@ int main (int argc, char ** argv)
// SU5::printGenerators();
// SU5::testGenerators();
Grid_finalize();
}

2
tests/Test_partfrac_force.cc
View File

@ -96,7 +96,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

2
tests/Test_rect_force.cc
View File

@ -81,7 +81,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

2
tests/Test_wilson_force.cc
View File

@ -93,7 +93,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
Hmom -= real(sum(trace(mommu*mommu)));

2
tests/Test_wilson_force_phiMdagMphi.cc
View File

@ -103,7 +103,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
// Dw.DoubleStore(Dw.Umu,Uprime); // update U _and_ Udag
Dw.DhopDirDisp(phi,Ftmp,mu,mu+4,DaggerYes);

2
tests/Test_wilson_force_phiMphi.cc
View File

@ -86,7 +86,7 @@ int main (int argc, char ** argv)
LatticeColourMatrix Umu_save(&Grid);
LatticeColourMatrix dU (&Grid);
LatticeColourMatrix mom(&Grid);
SU3::GaussianLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg
// check mom is as i expect