From 0bd296dda42b7cd65becf643fa62b4b1ad0be53d Mon Sep 17 00:00:00 2001
From: Guido Cossu <guido.cossu@ed.ac.uk>
Date: Wed, 14 Dec 2016 03:15:09 +0000
Subject: [PATCH] Adding check of the Dag part in the benchmark

---
 benchmarks/Benchmark_dwf.cc                  | 95 +++++++++++++++-----
 lib/algorithms/iterative/ConjugateGradient.h |  1 -
 lib/simd/Grid_avx.h                          |  4 +-
 3 files changed, 76 insertions(+), 24 deletions(-)

diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc
index 10e4521b..76625973 100644
--- a/benchmarks/Benchmark_dwf.cc
+++ b/benchmarks/Benchmark_dwf.cc
@@ -163,6 +163,24 @@ int main (int argc, char ** argv)
     Dw.Report();
   }
 
+  if (1) {  // Naive wilson dag implementation
+    ref = zero;
+    for (int mu = 0; mu < Nd; mu++) {
+      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
+      tmp = U[mu] * Cshift(src, mu + 1, 1);
+      for (int i = 0; i < ref._odata.size(); i++) {
+        ref._odata[i] += tmp._odata[i] + Gamma(Gmu[mu]) * tmp._odata[i];
+      }
+
+      tmp = adj(U[mu]) * src;
+      tmp = Cshift(tmp, mu + 1, -1);
+      for (int i = 0; i < ref._odata.size(); i++) {
+        ref._odata[i] += tmp._odata[i] - Gamma(Gmu[mu]) * tmp._odata[i];
+      }
+    }
+    ref = -0.5 * ref;
+  }
+
   if (1)
   {
 
@@ -245,6 +263,33 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage<<" difference between normal and simd is "<<sum<<std::endl;
     assert (sum< 1.0e-4 );
 
+    // Check Dag
+    std::cout << GridLogMessage << "Compare WilsonFermion5D<DomainWallVec5dImplR>::Dhop to naive wilson implementation Dag to verify correctness" << std::endl;
+    sDw.Dhop(ssrc,sresult,1);
+    sum=0;
+    for(int x=0;x<latt4[0];x++){
+    for(int y=0;y<latt4[1];y++){
+    for(int z=0;z<latt4[2];z++){
+    for(int t=0;t<latt4[3];t++){
+    for(int s=0;s<Ls;s++){
+      std::vector<int> site({s,x,y,z,t});
+      SpinColourVector normal, simd;
+      peekSite(normal,ref,site);
+      peekSite(simd,sresult,site);
+      sum=sum+norm2(normal-simd);
+      if (norm2(normal-simd) > 1.0e-6 ) {
+        std::cout << "site "<<x<<","<<y<<","<<z<<","<<t<<","<<s<<" "<<norm2(normal-simd)<<std::endl;
+        std::cout << "site "<<x<<","<<y<<","<<z<<","<<t<<","<<s<<" normal "<<normal<<std::endl;
+        std::cout << "site "<<x<<","<<y<<","<<z<<","<<t<<","<<s<<" simd   "<<simd<<std::endl;
+      }
+    }}}}}
+    std::cout<<GridLogMessage<<" difference between normal and simd is "<<sum<<std::endl;
+    assert (sum< 1.0e-4 );
+
+
+
+
+
 
     if (1) {
 
@@ -301,42 +346,50 @@ int main (int argc, char ** argv)
       ssrc_e = ssrc_e - sr_e;
       RealD error = norm2(ssrc_e);
 
-      std::cout<<GridLogMessage << "sE norm diff   "<< norm2(ssrc_e)<< "  vec nrm"<<norm2(sr_e) <<std::endl;
+      std::cout<<GridLogMessage << "sE norm diff   "<< norm2(ssrc_e)<< "  vec nrm: "<<norm2(sr_e) <<std::endl;
       ssrc_o = ssrc_o - sr_o;
 
       error+= norm2(ssrc_o);
-      std::cout<<GridLogMessage << "sO norm diff   "<< norm2(ssrc_o)<< "  vec nrm"<<norm2(sr_o) <<std::endl;
+      std::cout<<GridLogMessage << "sO norm diff   "<< norm2(ssrc_o)<< "  vec nrm: "<<norm2(sr_o) <<std::endl;
       if(error>1.0e-4) { 
 	setCheckerboard(ssrc,ssrc_o);
 	setCheckerboard(ssrc,ssrc_e);
 	std::cout<< ssrc << std::endl;
       }
+
+      // Check the dag
+      std::cout << GridLogMessage << "Compare WilsonFermion5D<DomainWallVec5dImplR>::DhopEO to Dhop to verify correctness" << std::endl;
+      pickCheckerboard(Even,ssrc_e,ssrc);
+      pickCheckerboard(Odd,ssrc_o,ssrc);
+      sDw.DhopEO(ssrc_o,sr_e,DaggerYes);
+      sDw.DhopOE(ssrc_e,sr_o,DaggerYes);
+      sDw.Dhop  (ssrc  ,sresult,DaggerYes);
+
+      pickCheckerboard(Even,ssrc_e,sresult);
+      pickCheckerboard(Odd ,ssrc_o,sresult);
+      ssrc_e = ssrc_e - sr_e;
+      error = norm2(ssrc_e);
+
+      std::cout<<GridLogMessage << "sE norm diff   "<< norm2(ssrc_e)<< "  vec nrm: "<<norm2(sr_e) <<std::endl;
+      ssrc_o = ssrc_o - sr_o;
+
+      error+= norm2(ssrc_o);
+      std::cout<<GridLogMessage << "sO norm diff   "<< norm2(ssrc_o)<< "  vec nrm: "<<norm2(sr_o) <<std::endl;
+      if(error>1.0e-4) { 
+        setCheckerboard(ssrc,ssrc_o);
+        setCheckerboard(ssrc,ssrc_e);
+        std::cout<< ssrc << std::endl;
+      } 
+
+
     }
 
 
   }
 
-  if (1)
-  { // Naive wilson dag implementation
-    ref = zero;
-    for(int mu=0;mu<Nd;mu++){
 
-      //    ref =  src - Gamma(Gamma::GammaX)* src ; // 1+gamma_x
-      tmp = U[mu]*Cshift(src,mu+1,1);
-      for(int i=0;i<ref._odata.size();i++){
-  ref._odata[i]+= tmp._odata[i] + Gamma(Gmu[mu])*tmp._odata[i]; ;
-      }
-
-      tmp =adj(U[mu])*src;
-      tmp =Cshift(tmp,mu+1,-1);
-      for(int i=0;i<ref._odata.size();i++){
-  ref._odata[i]+= tmp._odata[i] - Gamma(Gmu[mu])*tmp._odata[i]; ;
-      }
-    }
-    ref = -0.5*ref;
-  }
   Dw.Dhop(src,result,1);
-  std::cout << GridLogMessage << "Compare to naive wilson implementation Dag to verify correctness" << std::endl;
+  std::cout << GridLogMessage << "Compare DomainWallFermionR::Dhop to naive wilson implementation Dag to verify correctness" << std::endl;
   std::cout<<GridLogMessage << "Called DwDag"<<std::endl;
   std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
   std::cout<<GridLogMessage << "norm ref    "<< norm2(ref)<<std::endl;
diff --git a/lib/algorithms/iterative/ConjugateGradient.h b/lib/algorithms/iterative/ConjugateGradient.h
index 6fb69540..3bf3a5a3 100644
--- a/lib/algorithms/iterative/ConjugateGradient.h
+++ b/lib/algorithms/iterative/ConjugateGradient.h
@@ -134,7 +134,6 @@ class ConjugateGradient : public OperatorFunction<Field> {
       std::cout << GridLogDebug << "a = "<< a << " b_pred = "<< b_pred << "  b = "<< b << std::endl;
       std::cout << GridLogDebug << "qq = "<< qq << " d = "<< d << "  c = "<< c << std::endl;
 
-
       // Stopping condition
       if (cp <= rsq) {
         SolverTimer.Stop();
diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h
index 36360102..77927b10 100644
--- a/lib/simd/Grid_avx.h
+++ b/lib/simd/Grid_avx.h
@@ -514,7 +514,7 @@ namespace Optimization {
     template<int n>
     static inline __m256 tRotate(__m256 in){ 
       __m256 tmp = Permute::Permute0(in);
-      __m256 ret;
+      __m256 ret = in;
       if ( n > 3 ) { 
 	_mm256_alignr_epi32_grid(ret,in,tmp,n);  
       } else {
@@ -526,7 +526,7 @@ namespace Optimization {
     template<int n>
     static inline __m256d tRotate(__m256d in){ 
       __m256d tmp = Permute::Permute0(in);
-      __m256d ret;
+      __m256d ret = in;
       if ( n > 1 ) {
 	_mm256_alignr_epi64_grid(ret,in,tmp,n);          
       } else {