diff --git a/Grid/tensors/Tensor_exp.h b/Grid/tensors/Tensor_exp.h index 1f637d5f..45ae77fb 100644 --- a/Grid/tensors/Tensor_exp.h +++ b/Grid/tensors/Tensor_exp.h @@ -52,12 +52,17 @@ template accelerator_inline iVector Exponentiate(c return ret; } - - // Specialisation: Cayley-Hamilton exponential for SU(3) -#ifndef GRID_CUDA +#if 0 template::TensorLevel == 0>::type * =nullptr> -accelerator_inline iMatrix Exponentiate(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) +accelerator_inline iMatrix Exponentiated(const iMatrix &arg, RealD alpha , Integer Nexp = DEFAULT_MAT_EXP ) +{ + return ExponentiateCayleyHamilton(arg,alpha); +} +#endif + +template::TensorLevel == 0>::type * =nullptr> +accelerator_inline iMatrix ExponentiateCayleyHamilton(const iMatrix &arg, RealD alpha ) { // for SU(3) 2x faster than the std implementation using Nexp=12 // notice that it actually computes @@ -115,8 +120,6 @@ accelerator_inline iMatrix Exponentiate(const iMatrix &arg, Re return (f0 * unit + timesMinusI(f1) * arg*alpha - f2 * iQ2); } -#endif - // General exponential template::TensorLevel == 0 >::type * =nullptr> @@ -129,8 +132,8 @@ accelerator_inline iMatrix Exponentiate(const iMatrix &arg, Re typedef iMatrix mat; mat unit(1.0); mat temp(unit); - for(int i=Nexp; i>=1;--i){ - temp *= alpha/RealD(i); + for(int n=Nexp; n>=1;--n){ + temp *= alpha/RealD(n); temp = unit + temp*arg; } return temp;