portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 09:15:38 +01:00
Code Releases Activity

Using Cayley-Hamilton form for the exponential of SU(3) matrices

Browse Source
This commit is contained in:
Guido Cossu 2017-05-25 12:07:47 +01:00
parent 15e801af3f
commit ab3596d4d3
4 changed files with 106 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
lib/lattice/Lattice_unary.h
View File

@ -62,17 +62,12 @@ namespace Grid {
return ret;
}
template<class obj> Lattice<obj> expMat(const Lattice<obj> &rhs, ComplexD alpha, Integer Nexp = DEFAULT_MAT_EXP){
template<class obj> Lattice<obj> expMat(const Lattice<obj> &rhs, RealD alpha, Integer Nexp = DEFAULT_MAT_EXP){
Lattice<obj> ret(rhs._grid);
ret.checkerboard = rhs.checkerboard;
conformable(ret,rhs);
obj unit(1.0);
parallel_for(int ss=0;ss<rhs._grid->oSites();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;

15
lib/qcd/action/gauge/GaugeImplTypes.h
View File

@ -94,12 +94,17 @@ 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<LorentzIndex>(U, mu);
auto Pmu = PeekIndex<LorentzIndex>(P, mu);
Umu = expMat(Pmu, ep, Nexp) * Umu;
PokeIndex<LorentzIndex>(U, ProjectOnGroup(Umu), mu);
//static std::chrono::duration<double> diff;
//auto start = std::chrono::high_resolution_clock::now();
parallel_for(int ss=0;ss<P._grid->oSites();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){

2
lib/qcd/smearing/StoutSmearing.h
View File

@ -58,6 +58,8 @@ class Smear_Stout : public Smear<Gimpl> {
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

93
lib/tensors/Tensor_exp.h
View File

@ -37,30 +37,105 @@ namespace Grid {
///////////////////////////////////////////////
template<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP)
template<class vtype> inline iScalar<vtype> Exponentiate(const iScalar<vtype>&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP)
{
iScalar<vtype> ret;
ret._internal = Exponentiate(r._internal, alpha, Nexp);
return ret;
}
template<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &arg, ComplexD alpha , Integer Nexp = DEFAULT_MAT_EXP )
template<class vtype, int N> inline iVector<vtype, N> Exponentiate(const iVector<vtype,N>&r, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP)
{
iMatrix<vtype,N> unit(1.0);
iMatrix<vtype,N> temp(unit);
iVector<vtype, N> 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<class vtype, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0>::type * =nullptr>
inline iMatrix<vtype,3> Exponentiate(const iMatrix<vtype,3> &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<vtype,3> mat;
typedef iScalar<vtype> 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<class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
inline iMatrix<vtype,N> Exponentiate(const iMatrix<vtype,N> &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<vtype,N> mat;
mat unit(1.0);
mat temp(unit);
for(int i=Nexp; i>=1;--i){
temp *= alpha/ComplexD(i);
temp = unit + temp*arg;
temp *= alpha/ComplexD(i);
temp = unit + temp*arg;
}
return temp;
}
}
#endif