From ab3596d4d368c9112a3d169d3cd8459d0f072f05 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:07:47 +0100 Subject: [PATCH 1/4] Using Cayley-Hamilton form for the exponential of SU(3) matrices --- lib/lattice/Lattice_unary.h | 9 +-- lib/qcd/action/gauge/GaugeImplTypes.h | 21 ++++-- lib/qcd/smearing/StoutSmearing.h | 2 + lib/tensors/Tensor_exp.h | 103 ++++++++++++++++++++++---- 4 files changed, 106 insertions(+), 29 deletions(-) diff --git a/lib/lattice/Lattice_unary.h b/lib/lattice/Lattice_unary.h index 0afeaa32..44b7b4f1 100644 --- a/lib/lattice/Lattice_unary.h +++ b/lib/lattice/Lattice_unary.h @@ -62,17 +62,12 @@ namespace Grid { return ret; } - template Lattice expMat(const Lattice &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){ + template Lattice expMat(const Lattice &rhs, RealD alpha, Integer Nexp = DEFAULT_MAT_EXP){ Lattice ret(rhs._grid); ret.checkerboard = rhs.checkerboard; conformable(ret,rhs); - obj unit(1.0); parallel_for(int ss=0;ssoSites();ss++){ - //ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); - ret._odata[ss] = unit; - for(int i=Nexp; i>=1;--i) - ret._odata[ss] = unit + ret._odata[ss]*rhs._odata[ss]*(alpha/RealD(i)); - + ret._odata[ss]=Exponentiate(rhs._odata[ss],alpha, Nexp); } return ret; diff --git a/lib/qcd/action/gauge/GaugeImplTypes.h b/lib/qcd/action/gauge/GaugeImplTypes.h index 9d36eead..29e79581 100644 --- a/lib/qcd/action/gauge/GaugeImplTypes.h +++ b/lib/qcd/action/gauge/GaugeImplTypes.h @@ -59,7 +59,7 @@ public: typedef iImplGaugeLink SiteLink; typedef iImplGaugeField SiteField; - typedef Lattice LinkField; + typedef Lattice LinkField; typedef Lattice Field; // Guido: we can probably separate the types from the HMC functions @@ -80,7 +80,7 @@ public: /////////////////////////////////////////////////////////// // Move these to another class - // HMC auxiliary functions + // HMC auxiliary functions static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) { // specific for SU gauge fields LinkField Pmu(P._grid); @@ -92,14 +92,19 @@ public: } static inline Field projectForce(Field &P) { return Ta(P); } - + static inline void update_field(Field& P, Field& U, double ep){ - for (int mu = 0; mu < Nd; mu++) { - auto Umu = PeekIndex(U, mu); - auto Pmu = PeekIndex(P, mu); - Umu = expMat(Pmu, ep, Nexp) * Umu; - PokeIndex(U, ProjectOnGroup(Umu), mu); + //static std::chrono::duration diff; + + //auto start = std::chrono::high_resolution_clock::now(); + parallel_for(int ss=0;ssoSites();ss++){ + for (int mu = 0; mu < Nd; mu++) + U[ss]._internal[mu] = ProjectOnGroup(Exponentiate(P[ss]._internal[mu], ep, Nexp) * U[ss]._internal[mu]); } + + //auto end = std::chrono::high_resolution_clock::now(); + // diff += end - start; + // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n"; } static inline RealD FieldSquareNorm(Field& U){ diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h index b50ffe21..bfc37d0d 100644 --- a/lib/qcd/smearing/StoutSmearing.h +++ b/lib/qcd/smearing/StoutSmearing.h @@ -58,6 +58,8 @@ class Smear_Stout : public Smear { SmearBase->smear(C, U); }; + + // Repetion of code here (use the Tensor_exp.h function) void exponentiate_iQ(GaugeLinkField& e_iQ, const GaugeLinkField& iQ) const { // Put this outside // only valid for SU(3) matrices diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index ad0d2071..64c06bb5 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -37,30 +37,105 @@ namespace Grid { /////////////////////////////////////////////// - template inline iScalar Exponentiate(const iScalar&r, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP) + template inline iScalar Exponentiate(const iScalar&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP) { iScalar ret; ret._internal = Exponentiate(r._internal, alpha, Nexp); return ret; } - - template::TensorLevel == 0 >::type * =nullptr> - inline iMatrix Exponentiate(const iMatrix &arg, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP ) +template inline iVector Exponentiate(const iVector&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP) { - iMatrix unit(1.0); - iMatrix temp(unit); - - for(int i=Nexp; i>=1;--i){ - temp *= alpha/ComplexD(i); - temp = unit + temp*arg; - } - - return temp; - + iVector ret; + for (int i = 0; i < N; i++) + ret._internal[i] = Exponentiate(r._internal[i], alpha, Nexp); + return ret; } + // Specialisation: Cayley-Hamilton exponential for SU(3) + template::TensorLevel == 0>::type * =nullptr> + inline iMatrix Exponentiate(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) + { + // for SU(3) 2x faster than the std implementation using Nexp=12 + // notice that it actually computes + // exp ( input matrix ) + // the i sign is coming from outside + // input matrix is anti-hermitian NOT hermitian + typedef iMatrix mat; + typedef iScalar scalar; + mat unit(1.0); + mat temp(unit); + const Complex one_over_three = 1.0 / 3.0; + const Complex one_over_two = 1.0 / 2.0; + + scalar c0, c1, tmp, c0max, theta, u, w; + scalar xi0, u2, w2, cosw; + scalar fden, h0, h1, h2; + scalar e2iu, emiu, ixi0, qt; + scalar f0, f1, f2; + scalar unity(1.0); + + mat iQ2 = arg*arg*alpha*alpha; + mat iQ3 = arg*iQ2*alpha; + // sign in c0 from the conventions on the Ta + c0 = -imag( trace(iQ3) ) * one_over_three; + c1 = -real( trace(iQ2) ) * one_over_two; + + // Cayley Hamilton checks to machine precision, tested + tmp = c1 * one_over_three; + c0max = 2.0 * pow(tmp, 1.5); + + theta = acos(c0 / c0max) * one_over_three; + u = sqrt(tmp) * cos(theta); + w = sqrt(c1) * sin(theta); + + xi0 = sin(w) / w; + u2 = u * u; + w2 = w * w; + cosw = cos(w); + + ixi0 = timesI(xi0); + emiu = cos(u) - timesI(sin(u)); + e2iu = cos(2.0 * u) + timesI(sin(2.0 * u)); + + h0 = e2iu * (u2 - w2) + + emiu * ((8.0 * u2 * cosw) + (2.0 * u * (3.0 * u2 + w2) * ixi0)); + h1 = e2iu * (2.0 * u) - emiu * ((2.0 * u * cosw) - (3.0 * u2 - w2) * ixi0); + h2 = e2iu - emiu * (cosw + (3.0 * u) * ixi0); + + fden = unity / (9.0 * u2 - w2); // reals + f0 = h0 * fden; + f1 = h1 * fden; + f2 = h2 * fden; + + return (f0 * unit + timesMinusI(f1) * arg*alpha - f2 * iQ2); + } + + + +// General exponential +template::TensorLevel == 0 >::type * =nullptr> + inline iMatrix Exponentiate(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) + { + // notice that it actually computes + // exp ( input matrix ) + // the i sign is coming from outside + // input matrix is anti-hermitian NOT hermitian + typedef iMatrix mat; + mat unit(1.0); + mat temp(unit); + for(int i=Nexp; i>=1;--i){ + temp *= alpha/ComplexD(i); + temp = unit + temp*arg; + } + return temp; + + } + + + + } #endif From 3c112a7a25eebf5d1ae611f900d0d65057336fb3 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:09:00 +0100 Subject: [PATCH 2/4] Small correction to the general exp definition --- lib/tensors/Tensor_exp.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index 64c06bb5..2f584134 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -127,7 +127,7 @@ template::Tens mat unit(1.0); mat temp(unit); for(int i=Nexp; i>=1;--i){ - temp *= alpha/ComplexD(i); + temp *= alpha/ReadD(i); temp = unit + temp*arg; } return temp; From 75856f29454eaa3cfd7f002fa90ebd074ba663f4 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Thu, 25 May 2017 12:44:56 +0100 Subject: [PATCH 3/4] Compilation fix in the Tensor_exp --- lib/tensors/Tensor_exp.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/tensors/Tensor_exp.h b/lib/tensors/Tensor_exp.h index 2f584134..e18fed70 100644 --- a/lib/tensors/Tensor_exp.h +++ b/lib/tensors/Tensor_exp.h @@ -127,7 +127,7 @@ template::Tens mat unit(1.0); mat temp(unit); for(int i=Nexp; i>=1;--i){ - temp *= alpha/ReadD(i); + temp *= alpha/RealD(i); temp = unit + temp*arg; } return temp; From a74c34315c79c4b52871e87b312943174dbaf1e2 Mon Sep 17 00:00:00 2001 From: Antonin Portelli Date: Thu, 25 May 2017 14:27:49 +0100 Subject: [PATCH 4/4] Bootstrap script fix --- bootstrap.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bootstrap.sh b/bootstrap.sh index 66eb63ee..dfb6735d 100755 --- a/bootstrap.sh +++ b/bootstrap.sh @@ -1,4 +1,4 @@ -]#!/usr/bin/env bash +#!/usr/bin/env bash EIGEN_URL='http://bitbucket.org/eigen/eigen/get/3.3.3.tar.bz2'