From 8dd099267d657735bd45582d63f8a567806606c4 Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Sun, 3 Jul 2016 12:28:25 +0100
Subject: [PATCH] Corrected a bug in the Expression Templates (acso and asin
 were wrong)

---
 lib/lattice/Lattice_ET.h         |  4 +--
 lib/qcd/smearing/StoutSmearing.h | 61 +++++++++++++++++++++++++-------
 lib/tensors/Tensor_unary.h       |  2 ++
 3 files changed, 52 insertions(+), 15 deletions(-)

diff --git a/lib/lattice/Lattice_ET.h b/lib/lattice/Lattice_ET.h
index 77e0ce5d..98c38e96 100644
--- a/lib/lattice/Lattice_ET.h
+++ b/lib/lattice/Lattice_ET.h
@@ -219,8 +219,8 @@ inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3
  GridUnopClass(UnaryRsqrt,rsqrt(a));
  GridUnopClass(UnarySin,sin(a));
  GridUnopClass(UnaryCos,cos(a));
- GridUnopClass(UnaryAsin,sin(a));
- GridUnopClass(UnaryAcos,cos(a));
+ GridUnopClass(UnaryAsin,asin(a));
+ GridUnopClass(UnaryAcos,acos(a));
  GridUnopClass(UnaryLog,log(a));
  GridUnopClass(UnaryExp,exp(a));
  
diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h
index 6a061bac..3c323c47 100644
--- a/lib/qcd/smearing/StoutSmearing.h
+++ b/lib/qcd/smearing/StoutSmearing.h
@@ -30,6 +30,8 @@
   			~Smear_Stout(){}
 
   			void smear(GaugeField& u_smr,const GaugeField& U) const{
+
+
   				GaugeField C(U._grid);
   				GaugeLinkField tmp(U._grid), iq_mu(U._grid), Umu(U._grid);
 
@@ -41,9 +43,14 @@
   				{
   					tmp = peekLorentz(C,mu);
   					Umu = peekLorentz(U,mu);
+		  			std::cout << "source matrix " << Umu << std::endl;	
+
 		  			iq_mu = Ta(tmp * adj(Umu)); // iq_mu = Ta(Omega_mu) to match the signs with the paper
+
 		  			exponentiate_iQ(tmp, iq_mu);  
+		  			//Debug check
 		  			GaugeLinkField check = adj(tmp) * tmp - 1.0;
+		  			std::cout << "check " << check << std::endl;
 					pokeLorentz(u_smr, tmp*Umu, mu);// u_smr = exp(iQ_mu)*U_mu
 		  		}
 
@@ -90,9 +97,16 @@
 
 		  		set_uw_complex(u, w, iQ2, iQ3);
 		  		set_fj_complex(f0, f1, f2, u, w);
-		  		
+
+		  		std::cout << "f0 " << f0 << std::endl;
+		  		std::cout << "f1 " << f1 << std::endl;
+		  		std::cout << "f2 " << f2 << std::endl;
+	  			std::cout << "iQ " << iQ << std::endl;	
+	  			std::cout << "iQ2 " << iQ2 << std::endl;	
 
 		  		e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2;
+
+
 		  	};
 
 
@@ -128,13 +142,27 @@
 				// sign in c0 from the conventions on the Ta
 				c0    = - real(timesMinusI(trace(iQ3))) * one_over_three; //temporary hack 
 				c1    = - real(trace(iQ2)) * one_over_two;
+
+				//Cayley Hamilton checks to machine precision, tested
+
+				std::cout << "c0 " << c0 << std::endl;
+				std::cout << "c1 " << c1 << std::endl;
+
 				tmp   = c1 * one_over_three;
 				c0max = 2.0 * pow(tmp, 1.5);
 
-				theta = acos(c0/c0max);
+				std::cout << "c0max " << c0max << std::endl;
+				LatticeComplex tempratio = c0/c0max;
+				std::cout << "ratio c0/c0max " << tempratio << std::endl;
+				theta = acos(c0/c0max); // divide by three here, now leave as it is
+				std::cout << "theta " << theta << std::endl;
+
 				u = sqrt(tmp) * cos( theta * one_over_three);
 				w = sqrt(c1)  * sin( theta * one_over_three);
 
+				std::cout << "u " << u << std::endl;
+				std::cout << "w " << w << std::endl;
+
 			}
 
 
@@ -177,35 +205,42 @@
 				const LatticeComplex& u, const LatticeComplex& w) const{
 
 				GridBase *grid = u._grid;
-				LatticeComplex xi0(grid), u2(grid), w2(grid), cosw(grid), tmp(grid);
+				LatticeComplex xi0(grid), u2(grid), w2(grid), cosw(grid);
 				LatticeComplex fden(grid);
 				LatticeComplex h0(grid), h1(grid), h2(grid);
 				LatticeComplex e2iu(grid), emiu(grid), ixi0(grid), qt(grid);
+				LatticeComplex unity(grid);
+				unity = 1.0;
 
 				xi0   = sin(w)/w;//func_xi0(w);
+				std::cout << "xi0 " << xi0 << std::endl;
 				u2    = u * u;
+				std::cout << "u2 " << u2 << std::endl;
 				w2    = w * w;
+				std::cout << "w2 " << w2 << std::endl;
 				cosw  = cos(w);
+				std::cout << "cosw " << cosw << std::endl;
 
 				ixi0  = timesI(xi0);
-				emiu  = cos(u) - timesI(sin(u));
+				emiu  = cos(u)     - timesI(sin(u));
 				e2iu  = cos(2.0*u) + timesI(sin(2.0*u));
+				std::cout << "emiu " << emiu << std::endl;
+				std::cout << "e2iu " << e2iu << std::endl;
 
-				h0    = e2iu * (u2 - w2) + emiu *( (8.0*u2*cosw) +
-					(2.0*u*(3.0*u2 + w2)*ixi0));
+				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);
 
-				h1    = (2.0*u) * e2iu - emiu*( (2.0*u*cosw) -
-					(3.0*u2-w2)*ixi0);
+				std::cout << "h0 " << h0 << std::endl;
+				std::cout << "h1 " << h1 << std::endl;
+				std::cout << "h2 " << h2 << std::endl;
 
-				h2    = e2iu - emiu * (cosw + (3.0*u)*ixi0);
-
-				tmp   = 9.0*u2 - w2;
-				fden  = pow(tmp, -1.0);
+				fden   = unity/(9.0*u2 - w2);// reals
+				std::cout << "fden " << fden << std::endl;
 				f0    = h0 * fden;
 				f1    = h1 * fden;
 				f2    = h2 * fden;	
 
-
 			}
 
 
diff --git a/lib/tensors/Tensor_unary.h b/lib/tensors/Tensor_unary.h
index 1e1f09c0..dd05a4a7 100644
--- a/lib/tensors/Tensor_unary.h
+++ b/lib/tensors/Tensor_unary.h
@@ -86,6 +86,8 @@ UNARY(sqrt);
 UNARY(rsqrt);
 UNARY(sin);
 UNARY(cos);
+UNARY(asin);
+UNARY(acos);
 UNARY(log);
 UNARY(exp);
 UNARY(abs);