diff --git a/benchmarks/Benchmark_zmm.cc b/benchmarks/Benchmark_zmm.cc
index 8b2c548b..74674cbf 100644
--- a/benchmarks/Benchmark_zmm.cc
+++ b/benchmarks/Benchmark_zmm.cc
@@ -41,7 +41,7 @@ int main(int argc,char **argv)
   std::ofstream os("zmm.dat");
 
   os << "#V Ls Lxy Lzt C++ Asm OMP L1 " <<std::endl;
-  for(int L=8;L<=32;L+=4){
+  for(int L=4;L<=32;L+=4){
     for(int m=1;m<=2;m++){
       for(int Ls=8;Ls<=16;Ls+=8){
 	std::vector<int> grid({L,L,m*L,m*L});
diff --git a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
index 1e89fb00..c53d9318 100644
--- a/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
+++ b/lib/algorithms/iterative/ImplicitlyRestartedLanczos.h
@@ -564,6 +564,7 @@ until convergence
 	  
 	  for(int j=k1-1; j<k2+1; ++j){
 	    for(int k=0; k<Nm; ++k){
+	    B[j].checkerboard = evec[k].checkerboard;
 	      B[j] += Qt[k+Nm*j] * evec[k];
 	    }
 	  }
@@ -592,6 +593,7 @@ until convergence
 	  
 	  for(int j = 0; j<Nk; ++j){
 	    for(int k = 0; k<Nk; ++k){
+	    B[j].checkerboard = evec[k].checkerboard;
 	      B[j] += Qt[k+j*Nm] * evec[k];
 	    }
 //	    std::cout << "norm(B["<<j<<"])="<<norm2(B[j])<<std::endl;
diff --git a/lib/simd/Grid_avx.h b/lib/simd/Grid_avx.h
index f7292add..d05551c7 100644
--- a/lib/simd/Grid_avx.h
+++ b/lib/simd/Grid_avx.h
@@ -410,22 +410,22 @@ namespace Optimization {
   struct Permute{
 
     static inline __m256 Permute0(__m256 in){
-      return _mm256_permute2f128_ps(in,in,0x01);
+      return _mm256_permute2f128_ps(in,in,0x01); //ABCD EFGH -> EFGH ABCD
     };
     static inline __m256 Permute1(__m256 in){
-      return _mm256_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(1,0,3,2));
+      return _mm256_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(1,0,3,2)); //ABCD EFGH -> CDAB GHEF
     };
     static inline __m256 Permute2(__m256 in){
-      return _mm256_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(2,3,0,1));
+      return _mm256_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(2,3,0,1)); //ABCD EFGH -> BADC FEHG
     };
     static inline __m256 Permute3(__m256 in){
       return in;
     };
 
     static inline __m256d Permute0(__m256d in){
-      return _mm256_permute2f128_pd(in,in,0x01);
+      return _mm256_permute2f128_pd(in,in,0x01); //AB CD -> CD AB
     };
-    static inline __m256d Permute1(__m256d in){
+    static inline __m256d Permute1(__m256d in){ //AB CD -> BA DC
       return _mm256_shuffle_pd(in,in,0x5);
     };
     static inline __m256d Permute2(__m256d in){
diff --git a/lib/simd/Grid_empty.h b/lib/simd/Grid_empty.h
index 8858624e..5ab75de7 100644
--- a/lib/simd/Grid_empty.h
+++ b/lib/simd/Grid_empty.h
@@ -55,51 +55,67 @@ namespace Optimization {
   
   struct Vsplat{
     //Complex float
-    inline float operator()(float a, float b){
-      return 0;
+    inline u128f operator()(float a, float b){
+      u128f out; 
+      out.f[0] = a;
+      out.f[1] = b;
+      out.f[2] = a;
+      out.f[3] = b;
+      return out;
     }
     // Real float
-    inline float operator()(float a){
-      return 0;
+    inline u128f operator()(float a){
+      u128f out; 
+      out.f[0] = a;
+      out.f[1] = a;
+      out.f[2] = a;
+      out.f[3] = a;
+      return out;
     }
     //Complex double
-    inline double operator()(double a, double b){
-      return 0;
+    inline u128d operator()(double a, double b){
+      u128d out; 
+      out.f[0] = a;
+      out.f[1] = b;
+      return out;
     }
     //Real double
-    inline double operator()(double a){
-      return 0;
+    inline u128d operator()(double a){
+      u128d out; 
+      out.f[0] = a;
+      out.f[1] = a;
+      return out;
     }
     //Integer
     inline int operator()(Integer a){
-      return 0;
+      return a;
     }
   };
 
   struct Vstore{
     //Float 
-    inline void operator()(float a, float* F){
-      
+    inline void operator()(u128f a, float* F){
+      memcpy(F,a.f,4*sizeof(float));
     }
     //Double
-    inline void operator()(double a, double* D){
-     
+    inline void operator()(u128d a, double* D){
+      memcpy(D,a.f,2*sizeof(double));
     }
     //Integer
     inline void operator()(int a, Integer* I){
-      
+      I[0] = a;
     }
 
   };
 
   struct Vstream{
     //Float
-    inline void operator()(float * a, float b){
-     
+    inline void operator()(float * a, u128f b){
+      memcpy(a,b.f,4*sizeof(float));
     }
     //Double
-    inline void operator()(double * a, double b){
-     
+    inline void operator()(double * a, u128d b){
+      memcpy(a,b.f,2*sizeof(double));
     }
 
 
@@ -107,24 +123,40 @@ namespace Optimization {
 
   struct Vset{
     // Complex float 
-    inline float operator()(Grid::ComplexF *a){
-      return 0;
+    inline u128f operator()(Grid::ComplexF *a){
+      u128f out; 
+      out.f[0] = a[0].real();
+      out.f[1] = a[0].imag();
+      out.f[2] = a[1].real();
+      out.f[3] = a[1].imag();
+      return out;
     }
     // Complex double 
-    inline double operator()(Grid::ComplexD *a){
-      return 0;
+    inline u128d operator()(Grid::ComplexD *a){
+      u128d out; 
+      out.f[0] = a[0].real();
+      out.f[1] = a[0].imag();
+      return out;
     }
     // Real float 
-    inline float operator()(float *a){
-      return  0;
+    inline u128f operator()(float *a){
+      u128f out; 
+      out.f[0] = a[0];
+      out.f[1] = a[1];
+      out.f[2] = a[2];
+      out.f[3] = a[3];
+      return out;
     }
     // Real double
-    inline double operator()(double *a){
-      return 0;
+    inline u128d operator()(double *a){
+      u128d out; 
+      out.f[0] = a[0];
+      out.f[1] = a[1];
+      return out;
     }
     // Integer
     inline int operator()(Integer *a){
-      return 0;
+      return a[0];
     }
 
 
@@ -146,129 +178,198 @@ namespace Optimization {
   /////////////////////////////////////////////////////
   struct Sum{
     //Complex/Real float
-    inline float operator()(float a, float b){
-      return 0;
+    inline u128f operator()(u128f a, u128f b){
+      u128f out;
+      out.f[0] = a.f[0] + b.f[0];
+      out.f[1] = a.f[1] + b.f[1];
+      out.f[2] = a.f[2] + b.f[2];
+      out.f[3] = a.f[3] + b.f[3];
+      return out;
     }
     //Complex/Real double
-    inline double operator()(double a, double b){
-      return 0;
+    inline u128d operator()(u128d a, u128d b){
+      u128d out;
+      out.f[0] = a.f[0] + b.f[0];
+      out.f[1] = a.f[1] + b.f[1];
+      return out;
     }
     //Integer
     inline int operator()(int a, int b){
-      return 0;
+      return a + b;
     }
   };
 
   struct Sub{
     //Complex/Real float
-    inline float operator()(float a, float b){
-      return 0;
+    inline u128f operator()(u128f a, u128f b){
+      u128f out;
+      out.f[0] = a.f[0] - b.f[0];
+      out.f[1] = a.f[1] - b.f[1];
+      out.f[2] = a.f[2] - b.f[2];
+      out.f[3] = a.f[3] - b.f[3];
+      return out;
     }
     //Complex/Real double
-    inline double operator()(double a, double b){
-      return 0;
+    inline u128d operator()(u128d a, u128d b){
+      u128d out;
+      out.f[0] = a.f[0] - b.f[0];
+      out.f[1] = a.f[1] - b.f[1];
+      return out;
     }
     //Integer
     inline int operator()(int a, int b){
-      return 0;
+      return a-b;
     }
   };
 
   struct MultComplex{
     // Complex float
-    inline float operator()(float a, float b){
-      return 0;
+    inline u128f operator()(u128f a, u128f b){
+      u128f out;
+      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
+      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
+      out.f[2] = a.f[2]*b.f[2] - a.f[3]*b.f[3];
+      out.f[3] = a.f[2]*b.f[3] + a.f[3]*b.f[2];
+      return out;
     }
     // Complex double
-    inline double operator()(double a, double b){
-      return 0;
+    inline u128d operator()(u128d a, u128d b){
+      u128d out;
+      out.f[0] = a.f[0]*b.f[0] - a.f[1]*b.f[1];
+      out.f[1] = a.f[0]*b.f[1] + a.f[1]*b.f[0];
+      return out;
     }
   };
 
   struct Mult{
-    inline float  mac(float a, float b,double c){
-      return 0;
-    }
-    inline double mac(double a, double b,double c){
-      return 0;
-    }
+    //CK: Appear unneeded
+    // inline float  mac(float a, float b,double c){
+    //   return 0;
+    // }
+    // inline double mac(double a, double b,double c){
+    //   return 0;
+    // }
+
     // Real float
-    inline float operator()(float a, float b){
-      return 0;
+    inline u128f operator()(u128f a, u128f b){
+      u128f out;
+      out.f[0] = a.f[0]*b.f[0];
+      out.f[1] = a.f[1]*b.f[1];
+      out.f[2] = a.f[2]*b.f[2];
+      out.f[3] = a.f[3]*b.f[3];
+      return out;
     }
     // Real double
-    inline double operator()(double a, double b){
-      return 0;
+    inline u128d operator()(u128d a, u128d b){
+      u128d out;
+      out.f[0] = a.f[0]*b.f[0];
+      out.f[1] = a.f[1]*b.f[1];
+      return out;
     }
     // Integer
     inline int operator()(int a, int b){
-      return 0;
+      return a*b;
     }
   };
 
   struct Conj{
     // Complex single
-    inline float operator()(float in){
-      return 0;
+    inline u128f operator()(u128f in){
+      u128f out;
+      out.f[0] = in.f[0];
+      out.f[1] = -in.f[1];
+      out.f[2] = in.f[2];
+      out.f[3] = -in.f[3];
+      return out;
     }
     // Complex double
-    inline double operator()(double in){
-      return 0;
+    inline u128d operator()(u128d in){
+      u128d out;
+      out.f[0] = in.f[0];
+      out.f[1] = -in.f[1];
+      return out;
     }
     // do not define for integer input
   };
 
   struct TimesMinusI{
     //Complex single
-    inline float operator()(float in, float ret){
-      return 0;
+    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
+      u128f out;
+      out.f[0] = in.f[1];
+      out.f[1] = -in.f[0];
+      out.f[2] = in.f[3];
+      out.f[3] = -in.f[2];
+      return out;
     }
     //Complex double
-    inline double operator()(double in, double ret){
-      return 0;
+    inline u128d operator()(u128d in, u128d ret){
+      u128d out;
+      out.f[0] = in.f[1];
+      out.f[1] = -in.f[0];
+      return out;
     }
-
-
   };
 
   struct TimesI{
     //Complex single
-    inline float operator()(float in, float ret){
-      return 0;
+    inline u128f operator()(u128f in, u128f ret){ //note ret is ignored
+      u128f out;
+      out.f[0] = -in.f[1];
+      out.f[1] = in.f[0];
+      out.f[2] = -in.f[3];
+      out.f[3] = in.f[2];
+      return out;
     }
     //Complex double
-    inline double operator()(double in, double ret){
-      return 0;
+    inline u128d operator()(u128d in, u128d ret){
+      u128d out;
+      out.f[0] = -in.f[1];
+      out.f[1] = in.f[0];
+      return out;
     }
   };
 
   //////////////////////////////////////////////
   // Some Template specialization
   struct Permute{
-
-    static inline float Permute0(float in){
+    //We just have to mirror the permutes of Grid_sse4.h
+    static inline u128f Permute0(u128f in){ //AB CD -> CD AB
+      u128f out;
+      out.f[0] = in.f[2];
+      out.f[1] = in.f[3];
+      out.f[2] = in.f[0];
+      out.f[3] = in.f[1];
+      return out;
+    };
+    static inline u128f Permute1(u128f in){ //AB CD -> BA DC
+      u128f out;
+      out.f[0] = in.f[1];
+      out.f[1] = in.f[0];
+      out.f[2] = in.f[3];
+      out.f[3] = in.f[2];
+      return out;
+    };
+    static inline u128f Permute2(u128f in){
       return in;
     };
-    static inline float Permute1(float in){
-      return in;
-    };
-    static inline float Permute2(float in){
-      return in;
-    };
-    static inline float Permute3(float in){
+    static inline u128f Permute3(u128f in){
       return in;
     };
 
-    static inline double Permute0(double in){
+    static inline u128d Permute0(u128d in){ //AB -> BA
+      u128d out;
+      out.f[0] = in.f[1];
+      out.f[1] = in.f[0];
+      return out;      
+    };
+    static inline u128d Permute1(u128d in){
       return in;
     };
-    static inline double Permute1(double in){
+    static inline u128d Permute2(u128d in){
       return in;
     };
-    static inline double Permute2(double in){
-      return in;
-    };
-    static inline double Permute3(double in){
+    static inline u128d Permute3(u128d in){
       return in;
     };
 
@@ -280,26 +381,26 @@ namespace Optimization {
 
   //Complex float Reduce
   template<>
-  inline Grid::ComplexF Reduce<Grid::ComplexF, float>::operator()(float in){
-    return 0;
+  inline Grid::ComplexF Reduce<Grid::ComplexF, u128f>::operator()(u128f in){ //2 complex
+    return Grid::ComplexF(in.f[0] + in.f[2], in.f[1] + in.f[3]);
   }
   //Real float Reduce
   template<>
-  inline Grid::RealF Reduce<Grid::RealF, float>::operator()(float in){
-    return 0;
+  inline Grid::RealF Reduce<Grid::RealF, u128f>::operator()(u128f in){ //4 floats
+    return in.f[0] + in.f[1] + in.f[2] + in.f[3];
   }
   
   
   //Complex double Reduce
   template<>
-  inline Grid::ComplexD Reduce<Grid::ComplexD, double>::operator()(double in){
-    return 0;
+  inline Grid::ComplexD Reduce<Grid::ComplexD, u128d>::operator()(u128d in){ //1 complex
+    return Grid::ComplexD(in.f[0],in.f[1]);
   }
   
   //Real double Reduce
   template<>
-  inline Grid::RealD Reduce<Grid::RealD, double>::operator()(double in){
-    return 0;
+  inline Grid::RealD Reduce<Grid::RealD, u128d>::operator()(u128d in){ //2 doubles
+    return in.f[0] + in.f[1];
   }
 
   //Integer Reduce
@@ -314,8 +415,8 @@ namespace Optimization {
 //////////////////////////////////////////////////////////////////////////////////////
 // Here assign types 
 
-  typedef float SIMD_Ftype;  // Single precision type
-  typedef double SIMD_Dtype; // Double precision type
+  typedef Optimization::u128f SIMD_Ftype;  // Single precision type
+  typedef Optimization::u128d SIMD_Dtype; // Double precision type
   typedef int SIMD_Itype; // Integer type
 
   // prefetch utilities
diff --git a/lib/simd/Grid_sse4.h b/lib/simd/Grid_sse4.h
index 8abd4581..2af91f44 100644
--- a/lib/simd/Grid_sse4.h
+++ b/lib/simd/Grid_sse4.h
@@ -267,10 +267,10 @@ namespace Optimization {
   struct Permute{
 
     static inline __m128 Permute0(__m128 in){
-      return _mm_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(1,0,3,2));
+      return _mm_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(1,0,3,2)); //AB CD -> CD AB
     };
     static inline __m128 Permute1(__m128 in){
-      return _mm_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(2,3,0,1));
+      return _mm_shuffle_ps(in,in,_MM_SELECT_FOUR_FOUR(2,3,0,1)); //AB CD -> BA DC
     };
     static inline __m128 Permute2(__m128 in){
       return in;
@@ -279,7 +279,7 @@ namespace Optimization {
       return in;
     };
 
-    static inline __m128d Permute0(__m128d in){
+    static inline __m128d Permute0(__m128d in){ //AB -> BA
       return _mm_shuffle_pd(in,in,0x1);
     };
     static inline __m128d Permute1(__m128d in){
diff --git a/tests/Make.inc b/tests/Make.inc
index f7c83671..98323b45 100644
--- a/tests/Make.inc
+++ b/tests/Make.inc
@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gp_rect_force Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
+bin_PROGRAMS += Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gp_rect_force Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd 
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 5dff7e25..bbc2613c 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -20,6 +20,8 @@ endif
 
 if BUILD_ZMM
   bin_PROGRAMS=Test_zmm
+else
+  bin_PROGRAMS=
 endif
 
 include Make.inc