From ae9e248c95be15afd376dfed52f0374c3ff6fca0 Mon Sep 17 00:00:00 2001
From: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Date: Tue, 26 Feb 2019 11:29:12 +0000
Subject: [PATCH 1/8] Smearing

---
 Grid/Grid.h | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Grid/Grid.h b/Grid/Grid.h
index 9dcc207b..70859c81 100644
--- a/Grid/Grid.h
+++ b/Grid/Grid.h
@@ -42,6 +42,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/GridQCDcore.h>
 #include <Grid/qcd/action/Action.h>
 #include <Grid/qcd/utils/GaugeFix.h>
+#include <Grid/qcd/utils/CovariantSmearing.h>
 #include <Grid/qcd/smearing/Smearing.h>
 #include <Grid/parallelIO/MetaData.h>
 #include <Grid/qcd/hmc/HMC_aggregate.h>

From 8f661f6c05bc4ad7117d39f74e7ca8b7cf297ca8 Mon Sep 17 00:00:00 2001
From: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Date: Tue, 26 Feb 2019 11:31:00 +0000
Subject: [PATCH 2/8] Smearing for quark observables

---
 Grid/qcd/utils/CovariantSmearing.h | 87 ++++++++++++++++++++++++++++++
 1 file changed, 87 insertions(+)
 create mode 100644 Grid/qcd/utils/CovariantSmearing.h

diff --git a/Grid/qcd/utils/CovariantSmearing.h b/Grid/qcd/utils/CovariantSmearing.h
new file mode 100644
index 00000000..eae5ed71
--- /dev/null
+++ b/Grid/qcd/utils/CovariantSmearing.h
@@ -0,0 +1,87 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h
+
+Copyright (C) 2016
+
+Author: Azusa Yamaguchi
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template <class Gimpl> class CovariantSmearing : public Gimpl 
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  typedef typename Gimpl::GaugeLinkField GaugeMat;
+  typedef typename Gimpl::GaugeField GaugeLorentz;
+
+  template<typename T>
+  static void GaussianSmear(const std::vector<LatticeColourMatrix>& U, 
+			    T& chi, 
+			    const Real& width, int Iterations, int orthog)
+  {
+    GridBase *grid = chi._grid;
+    T psi(grid);
+
+    ////////////////////////////////////////////////////////////////////////////////////
+    // Follow Chroma conventions for width to keep compatibility with previous data
+    // Free field iterates 
+    //   chi = (1 - w^2/4N p^2)^N chi
+    //
+    //       ~ (e^(-w^2/4N p^2)^N chi
+    //       ~ (e^(-w^2/4 p^2) chi
+    //       ~ (e^(-w'^2/2 p^2) chi          [ w' = w/sqrt(2) ]
+    //
+    // Which in coordinate space is proportional to
+    //
+    //   e^(-x^2/w^2) = e^(-x^2/2w'^2) 
+    //
+    // The 4 is a bit unconventional from Gaussian width perspective, but... it's Chroma convention.
+    // 2nd derivative approx d^2/dx^2  =  x+mu + x-mu - 2x
+    //
+    // d^2/dx^2 = - p^2
+    //
+    // chi = ( 1 + w^2/4N d^2/dx^2 )^N chi
+    //
+    ////////////////////////////////////////////////////////////////////////////////////
+    Real coeff = (width*width) / Real(4*Iterations);
+  
+    int dims = Nd;
+    if( orthog < Nd ) dims=Nd-1;
+
+    for(int n = 0; n < Iterations; ++n) {
+      psi = (-2.0*dims)*chi;
+      for(int mu=0;mu<Nd;mu++) {
+	if ( mu != orthog ) { 
+	  psi = psi + Gimpl::CovShiftForward(U[mu],mu,chi);    
+	  psi = psi + Gimpl::CovShiftBackward(U[mu],mu,chi);    
+	}
+      }
+      chi = chi + coeff*psi;
+    }
+  }
+};
+}}

From 1ea64b24fe83bedb19219ea4368712f355bc5555 Mon Sep 17 00:00:00 2001
From: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Date: Tue, 26 Feb 2019 11:31:17 +0000
Subject: [PATCH 3/8] Smearing test. Test on free field.

---
 tests/smearing/Test_smearing.cc | 73 +++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 tests/smearing/Test_smearing.cc

diff --git a/tests/smearing/Test_smearing.cc b/tests/smearing/Test_smearing.cc
new file mode 100644
index 00000000..996671f0
--- /dev/null
+++ b/tests/smearing/Test_smearing.cc
@@ -0,0 +1,73 @@
+/*
+ *
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/Test_wilson_cg_prec.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          pRNG(&Grid);  pRNG.SeedFixedIntegers(seeds);
+
+  LatticeFermion src(&Grid); random(pRNG,src);
+  RealD nrm = norm2(src);
+  LatticeFermion result(&Grid); result=zero;
+  LatticeGaugeField Umu(&Grid);
+  //  SU3::HotConfiguration(pRNG,Umu);
+  SU3::ColdConfiguration(Umu);
+  std::vector<LatticeColourMatrix> U(4,&Grid);
+
+  for(int mu=0;mu<Nd;mu++){
+    U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+  }
+
+  std::vector<int> site({4,4,0,0});
+  src=zero;
+  SpinColourVector scv;
+  scv=zero;
+  scv()(0)(0) = 1.0;
+  pokeSite(scv,src,site);
+
+  CovariantSmearing<PeriodicGimplR>::GaussianSmear(U, src, 2.0, 50, Tdir);
+
+  std::cout << src <<std::endl;
+  
+}
+

From bdf87bc994de86ec354d44adbb5ff9f5d5f54e24 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 27 Feb 2019 02:27:09 +0000
Subject: [PATCH 4/8] Hadrons: beware of the nasty uninitialised twists

---
 Hadrons/Modules/MAction/DWF.hpp          | 24 ++++++++++++++++++----
 Hadrons/Modules/MAction/MobiusDWF.hpp    | 26 +++++++++++++++++++-----
 Hadrons/Modules/MAction/ScaledDWF.hpp    | 24 ++++++++++++++++++----
 Hadrons/Modules/MAction/Wilson.hpp       | 24 ++++++++++++++++++----
 Hadrons/Modules/MAction/WilsonClover.hpp | 25 +++++++++++++++++++----
 Hadrons/Modules/MAction/ZMobiusDWF.hpp   | 25 ++++++++++++++++++-----
 6 files changed, 122 insertions(+), 26 deletions(-)

diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp
index bac70ae8..21024a4f 100644
--- a/Hadrons/Modules/MAction/DWF.hpp
+++ b/Hadrons/Modules/MAction/DWF.hpp
@@ -112,8 +112,6 @@ void TDWF<FImpl>::setup(void)
                  << par().mass << ", M5= " << par().M5 << " and Ls= "
                  << par().Ls << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -121,8 +119,26 @@ void TDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename DomainWallFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
 }
diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp
index 6b960e2e..3ce02b91 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/MobiusDWF.hpp
@@ -112,17 +112,33 @@ void TMobiusDWF<FImpl>::setup(void)
                  << ", b= " << par().b << ", c= " << par().c
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
-
+                 
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename MobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                      implParams);
diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp
index 8742a820..3e052225 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.hpp
+++ b/Hadrons/Modules/MAction/ScaledDWF.hpp
@@ -111,8 +111,6 @@ void TScaledDWF<FImpl>::setup(void)
                  << ", scale= " << par().scale
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -120,8 +118,26 @@ void TScaledDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename ScaledShamirFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().scale,
                      implParams);
diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp
index 3ce5569e..11df89c3 100644
--- a/Hadrons/Modules/MAction/Wilson.hpp
+++ b/Hadrons/Modules/MAction/Wilson.hpp
@@ -109,15 +109,31 @@ void TWilson<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
 }
diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp
index a9ea8de5..d80bccdf 100644
--- a/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -112,17 +112,34 @@ void TWilsonClover<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
     LOG(Message) << "Clover term csw_r: " << par().csw_r
                  << " csw_t: " << par().csw_t
                  << std::endl;
+                 
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonCloverFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                      gridRb, par().mass, par().csw_r, par().csw_t, 
                      par().clover_anisotropy, implParams); 
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
index e780bd73..56efe1a2 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
@@ -118,10 +118,7 @@ void TZMobiusDWF<FImpl>::setup(void)
     {
         LOG(Message) << "  omega[" << i << "]= " << par().omega[i] << std::endl;
     }
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
-    env().createGrid(par().Ls);
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
     auto &grb4 = *envGetRbGrid(FermionField);
@@ -129,8 +126,26 @@ void TZMobiusDWF<FImpl>::setup(void)
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
     typename ZMobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, omega,
                      par().b, par().c, implParams);

From 7852181c2c532ab25ab2ae69b1ad71ce5df56482 Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Wed, 27 Feb 2019 02:27:40 +0000
Subject: [PATCH 5/8] Hadrons: uninitialised pointer fix (might have been
 harmless)

---
 Hadrons/Environment.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Hadrons/Environment.hpp b/Hadrons/Environment.hpp
index 05d3d8dd..24f6b31e 100644
--- a/Hadrons/Environment.hpp
+++ b/Hadrons/Environment.hpp
@@ -182,7 +182,7 @@ private:
     std::map<CoarseGridKey, GridPt>     gridCoarse5d_;
     unsigned int                        nd_;
     // random number generator
-    RngPt                               rng4d_;
+    RngPt                               rng4d_{nullptr};
     // object store
     std::vector<ObjInfo>                object_;
     std::map<std::string, unsigned int> objectAddress_;

From b7db99967a1cbbe84964caea41606f3d523aa74f Mon Sep 17 00:00:00 2001
From: Michael Marshall <michael.marshall@ed.ac.uk>
Date: Thu, 28 Feb 2019 20:06:59 +0000
Subject: [PATCH 6/8] Recommendations for Traits classes

---
 Grid/simd/Grid_vector_types.h |  20 ++--
 Grid/tensors/Tensor_class.h   | 100 +++++++++++++++++-
 Grid/tensors/Tensor_traits.h  | 194 ++++++++++++++++------------------
 3 files changed, 201 insertions(+), 113 deletions(-)

diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h
index c67e74cb..b34e056d 100644
--- a/Grid/simd/Grid_vector_types.h
+++ b/Grid/simd/Grid_vector_types.h
@@ -89,17 +89,25 @@ template <typename Condition, typename ReturnType> using NotEnableIf = Invoke<st
 ////////////////////////////////////////////////////////
 // Check for complexity with type traits
 template <typename T> struct is_complex : public std::false_type {};
-template <> struct is_complex<std::complex<double> > : public std::true_type {};
-template <> struct is_complex<std::complex<float> > : public std::true_type {};
+template <> struct is_complex<ComplexD> : public std::true_type {};
+template <> struct is_complex<ComplexF> : public std::true_type {};
 
-template <typename T>              using IfReal    = Invoke<std::enable_if<std::is_floating_point<T>::value, int> >;
+template<typename T, typename V=void> struct is_real : public std::false_type {};
+template<typename T> struct is_real<T, typename std::enable_if<std::is_floating_point<T>::value,
+  void>::type> : public std::true_type {};
+
+template<typename T, typename V=void> struct is_integer : public std::false_type {};
+template<typename T> struct is_integer<T, typename std::enable_if<std::is_integral<T>::value,
+  void>::type> : public std::true_type {};
+  
+template <typename T>              using IfReal    = Invoke<std::enable_if<is_real<T>::value, int> >;
 template <typename T>              using IfComplex = Invoke<std::enable_if<is_complex<T>::value, int> >;
-template <typename T>              using IfInteger = Invoke<std::enable_if<std::is_integral<T>::value, int> >;
+template <typename T>              using IfInteger = Invoke<std::enable_if<is_integer<T>::value, int> >;
 template <typename T1,typename T2> using IfSame    = Invoke<std::enable_if<std::is_same<T1,T2>::value, int> >;
 
-template <typename T>              using IfNotReal    = Invoke<std::enable_if<!std::is_floating_point<T>::value, int> >;
+template <typename T>              using IfNotReal    = Invoke<std::enable_if<!is_real<T>::value, int> >;
 template <typename T>              using IfNotComplex = Invoke<std::enable_if<!is_complex<T>::value, int> >;
-template <typename T>              using IfNotInteger = Invoke<std::enable_if<!std::is_integral<T>::value, int> >;
+template <typename T>              using IfNotInteger = Invoke<std::enable_if<!is_integer<T>::value, int> >;
 template <typename T1,typename T2> using IfNotSame    = Invoke<std::enable_if<!std::is_same<T1,T2>::value, int> >;
 
 ////////////////////////////////////////////////////////
diff --git a/Grid/tensors/Tensor_class.h b/Grid/tensors/Tensor_class.h
index c7f868db..c87503c6 100644
--- a/Grid/tensors/Tensor_class.h
+++ b/Grid/tensors/Tensor_class.h
@@ -62,7 +62,7 @@ class iScalar {
   // get double precision version
   typedef iScalar<typename GridTypeMapper<vtype>::DoublePrecision> DoublePrecision;
   
-  enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
+  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
 
   // Scalar no action
   //  template<int Level> using tensor_reduce_level = typename
@@ -173,7 +173,37 @@ class iScalar {
     return stream;
   };
 
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return &_internal; }
 
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return _internal.begin(); }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return (&_internal) + 1; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return &_internal; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return _internal.begin(); }
+  
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return (&_internal) + 1; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
 };
 ///////////////////////////////////////////////////////////
 // Allows to turn scalar<scalar<scalar<double>>>> back to double.
@@ -218,7 +248,7 @@ class iVector {
     return *this;
   }
 
-  enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
+  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
   iVector(const Zero &z) { *this = zero; };
   iVector() = default;
   /*
@@ -303,6 +333,38 @@ class iVector {
   //    strong_inline vtype && operator ()(int i) {
   //      return _internal[i];
   //    }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return _internal; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return _internal[0].begin(); }
+  
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return _internal + N; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return _internal; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return _internal[0].begin(); }
+  
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return _internal + N; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
 };
 
 template <class vtype, int N>
@@ -328,7 +390,7 @@ class iMatrix {
   typedef iScalar<tensor_reduced_v> tensor_reduced;
   typedef iMatrix<recurse_scalar_object, N> scalar_object;
 
-  enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1 };
+  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
 
   iMatrix(const Zero &z) { *this = zero; };
   iMatrix() = default;
@@ -458,6 +520,38 @@ class iMatrix {
   //  strong_inline vtype && operator ()(int i,int j) {
   //    return _internal[i][j];
   //  }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return _internal[0]; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline begin() const { return _internal[0][0].begin(); }
+  
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return _internal[0] + N * N; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
+  strong_inline end() const { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
+
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return _internal[0]; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline begin() { return _internal[0][0].begin(); }
+  
+  template <typename T = vtype>
+  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return _internal[0] + N * N; }
+  
+  template <typename T = vtype>
+  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
+  strong_inline end() { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
 };
 
 template <class v>
diff --git a/Grid/tensors/Tensor_traits.h b/Grid/tensors/Tensor_traits.h
index c1ef397a..d4854768 100644
--- a/Grid/tensors/Tensor_traits.h
+++ b/Grid/tensors/Tensor_traits.h
@@ -26,6 +26,21 @@ Author: Christopher Kelly <ckelly@phys.columbia.edu>
 
 namespace Grid {
 
+  // Forward declarations
+  template<class T>        class iScalar;
+  template<class T, int N> class iVector;
+  template<class T, int N> class iMatrix;
+
+  // These are the Grid tensors
+  template<typename T> struct isGridTensor
+  : public std::false_type { static constexpr bool notvalue = true; };
+  template<class T> struct isGridTensor<iScalar<T>>
+  : public std::true_type { static constexpr bool notvalue = false; };
+  template<class T, int N> struct isGridTensor<iVector<T, N>>
+  : public std::true_type { static constexpr bool notvalue = false; };
+  template<class T, int N> struct isGridTensor<iMatrix<T, N>>
+  : public std::true_type { static constexpr bool notvalue = false; };
+
 //////////////////////////////////////////////////////////////////////////////////
 // Want to recurse: GridTypeMapper<Matrix<vComplexD> >::scalar_type == ComplexD.
 // Use of a helper class like this allows us to template specialise and "dress"
@@ -40,25 +55,31 @@ namespace Grid {
 // to study C++11's type_traits.h file. (std::enable_if<isGridTensorType<vtype> >)
 //
 //////////////////////////////////////////////////////////////////////////////////
-  
-  template <class T> class GridTypeMapper {
-  public:
-    typedef typename T::scalar_type scalar_type;
-    typedef typename T::vector_type vector_type;
-    typedef typename T::vector_typeD vector_typeD;
-    typedef typename T::tensor_reduced tensor_reduced;
-    typedef typename T::scalar_object scalar_object;
-    typedef typename T::Complexified Complexified;
-    typedef typename T::Realified Realified;
-    typedef typename T::DoublePrecision DoublePrecision;
-    enum { TensorLevel = T::TensorLevel };
+
+  // This saves repeating common properties for supported Grid Scalar types
+  template<typename T, typename V=void> struct GridTypeMapper_Base {};
+  // TensorLevel    How many nested grid tensors
+  // Rank           Rank of the grid tensor
+  // count          Total number of elements, i.e. product of dimensions
+  // scalar_size    Size of the underlying fundamental object (tensor_reduced) in bytes
+  // size           Total size of all elements in bytes
+  // Dimension(dim) Size of dimension dim
+  template<typename T> struct GridTypeMapper_Base<T> {
+    static constexpr int TensorLevel = 0;
+    static constexpr int Rank = 0;
+    static constexpr std::size_t count = 1;
+    static constexpr std::size_t scalar_size = sizeof(T);
+    static constexpr std::size_t size = scalar_size * count;
+    static constexpr int Dimension(unsigned int dim) { return 0; }
   };
 
 //////////////////////////////////////////////////////////////////////////////////
 // Recursion stops with these template specialisations
 //////////////////////////////////////////////////////////////////////////////////
-  template<> class GridTypeMapper<RealF> {
-  public:
+
+  template<typename T> struct GridTypeMapper {};
+
+  template<> struct GridTypeMapper<RealF> : public GridTypeMapper_Base<RealF> {
     typedef RealF scalar_type;
     typedef RealF vector_type;
     typedef RealD vector_typeD;
@@ -67,10 +88,8 @@ namespace Grid {
     typedef ComplexF Complexified;
     typedef RealF Realified;
     typedef RealD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<RealD> {
-  public:
+  template<> struct GridTypeMapper<RealD> : public GridTypeMapper_Base<RealD> {
     typedef RealD scalar_type;
     typedef RealD vector_type;
     typedef RealD vector_typeD;
@@ -79,10 +98,8 @@ namespace Grid {
     typedef ComplexD Complexified;
     typedef RealD Realified;
     typedef RealD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<ComplexF> {
-  public:
+  template<> struct GridTypeMapper<ComplexF> : public GridTypeMapper_Base<ComplexF> {
     typedef ComplexF scalar_type;
     typedef ComplexF vector_type;
     typedef ComplexD vector_typeD;
@@ -91,10 +108,8 @@ namespace Grid {
     typedef ComplexF Complexified;
     typedef RealF Realified;
     typedef ComplexD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<ComplexD> {
-  public:
+  template<> struct GridTypeMapper<ComplexD> : public GridTypeMapper_Base<ComplexD> {
     typedef ComplexD scalar_type;
     typedef ComplexD vector_type;
     typedef ComplexD vector_typeD;
@@ -103,10 +118,8 @@ namespace Grid {
     typedef ComplexD Complexified;
     typedef RealD Realified;
     typedef ComplexD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<Integer> {
-  public:
+  template<> struct GridTypeMapper<Integer> : public GridTypeMapper_Base<Integer> {
     typedef Integer scalar_type;
     typedef Integer vector_type;
     typedef Integer vector_typeD;
@@ -115,11 +128,9 @@ namespace Grid {
     typedef void Complexified;
     typedef void Realified;
     typedef void DoublePrecision;
-    enum { TensorLevel = 0 };
   };
 
-  template<> class GridTypeMapper<vRealF> {
-  public:
+  template<> struct GridTypeMapper<vRealF> : public GridTypeMapper_Base<vRealF> {
     typedef RealF  scalar_type;
     typedef vRealF vector_type;
     typedef vRealD vector_typeD;
@@ -128,10 +139,8 @@ namespace Grid {
     typedef vComplexF Complexified;
     typedef vRealF Realified;
     typedef vRealD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<vRealD> {
-  public:
+  template<> struct GridTypeMapper<vRealD> : public GridTypeMapper_Base<vRealD> {
     typedef RealD  scalar_type;
     typedef vRealD vector_type;
     typedef vRealD vector_typeD;
@@ -140,10 +149,8 @@ namespace Grid {
     typedef vComplexD Complexified;
     typedef vRealD Realified;
     typedef vRealD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<vComplexH> {
-  public:
+  template<> struct GridTypeMapper<vComplexH> : public GridTypeMapper_Base<vComplexH> {
     typedef ComplexF  scalar_type;
     typedef vComplexH vector_type;
     typedef vComplexD vector_typeD;
@@ -152,10 +159,8 @@ namespace Grid {
     typedef vComplexH Complexified;
     typedef vRealH Realified;
     typedef vComplexD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<vComplexF> {
-  public:
+  template<> struct GridTypeMapper<vComplexF> : public GridTypeMapper_Base<vComplexF> {
     typedef ComplexF  scalar_type;
     typedef vComplexF vector_type;
     typedef vComplexD vector_typeD;
@@ -164,10 +169,8 @@ namespace Grid {
     typedef vComplexF Complexified;
     typedef vRealF Realified;
     typedef vComplexD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<vComplexD> {
-  public:
+  template<> struct GridTypeMapper<vComplexD> : public GridTypeMapper_Base<vComplexD> {
     typedef ComplexD  scalar_type;
     typedef vComplexD vector_type;
     typedef vComplexD vector_typeD;
@@ -176,10 +179,8 @@ namespace Grid {
     typedef vComplexD Complexified;
     typedef vRealD Realified;
     typedef vComplexD DoublePrecision;
-    enum { TensorLevel = 0 };
   };
-  template<> class GridTypeMapper<vInteger> {
-  public:
+  template<> struct GridTypeMapper<vInteger> : public GridTypeMapper_Base<vInteger> {
     typedef  Integer scalar_type;
     typedef vInteger vector_type;
     typedef vInteger vector_typeD;
@@ -188,57 +189,49 @@ namespace Grid {
     typedef void Complexified;
     typedef void Realified;
     typedef void DoublePrecision;
-    enum { TensorLevel = 0 };
   };
 
-  // First some of my own traits
-  template<typename T> struct isGridTensor {
-    static const bool value = true;
-    static const bool notvalue = false;
+#define GridTypeMapper_RepeatedTypes \
+  typedef typename ObjectTraits::scalar_type scalar_type; \
+  typedef typename ObjectTraits::vector_type vector_type; \
+  typedef typename ObjectTraits::vector_typeD vector_typeD; \
+  typedef typename ObjectTraits::tensor_reduced tensor_reduced; \
+  typedef typename ObjectTraits::scalar_object scalar_object; \
+  typedef typename ObjectTraits::Complexified Complexified; \
+  typedef typename ObjectTraits::Realified Realified; \
+  typedef typename ObjectTraits::DoublePrecision DoublePrecision; \
+  static constexpr int TensorLevel = BaseTraits::TensorLevel + 1; \
+  static constexpr std::size_t scalar_size = BaseTraits::scalar_size; \
+  static constexpr std::size_t size = scalar_size * count
+
+  template<typename T> struct GridTypeMapper<iScalar<T>> {
+    using ObjectTraits = iScalar<T>;
+    using BaseTraits = GridTypeMapper<T>;
+    static constexpr int Rank = 1 + BaseTraits::Rank;
+    static constexpr std::size_t count = 1 * BaseTraits::count;
+    static constexpr int Dimension(unsigned int dim) {
+      return ( dim == 0 ) ? 1 : BaseTraits::Dimension(dim - 1); }
+    GridTypeMapper_RepeatedTypes;
   };
-  template<> struct isGridTensor<int > {
-    static const bool value = false;
-    static const bool notvalue = true;
+
+  template<typename T, int N> struct GridTypeMapper<iVector<T, N>> {
+    using ObjectTraits = iVector<T, N>;
+    using BaseTraits = GridTypeMapper<T>;
+    static constexpr int Rank = 1 + BaseTraits::Rank;
+    static constexpr std::size_t count = N * BaseTraits::count;
+    static constexpr int Dimension(unsigned int dim) {
+      return ( dim == 0 ) ? N : BaseTraits::Dimension(dim - 1); }
+    GridTypeMapper_RepeatedTypes;
   };
-  template<> struct isGridTensor<RealD > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<RealF > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<ComplexD > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<ComplexF > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<Integer > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<vRealD > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<vRealF > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<vComplexD > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<vComplexF > {
-    static const bool value = false;
-    static const bool notvalue = true;
-  };
-  template<> struct isGridTensor<vInteger > {
-    static const bool value = false;
-    static const bool notvalue = true;
+
+  template<typename T, int N> struct GridTypeMapper<iMatrix<T, N>> {
+    using ObjectTraits = iMatrix<T, N>;
+    using BaseTraits = GridTypeMapper<T>;
+    static constexpr int Rank = 2 + BaseTraits::Rank;
+    static constexpr std::size_t count = N * N * BaseTraits::count;
+    static constexpr int Dimension(unsigned int dim) {
+      return ( dim == 0 || dim == 1 ) ? N : BaseTraits::Dimension(dim - 2); }
+    GridTypeMapper_RepeatedTypes;
   };
 
   // Match the index
@@ -263,20 +256,13 @@ namespace Grid {
     typedef T type;
   };
   
-  //Query if a tensor or Lattice<Tensor> is SIMD vector or scalar
-  template<typename T>
-  class isSIMDvectorized{
-    template<typename U>
-    static typename std::enable_if< !std::is_same< typename GridTypeMapper<typename getVectorType<U>::type>::scalar_type,   
-      typename GridTypeMapper<typename getVectorType<U>::type>::vector_type>::value, char>::type test(void *);
+  //Query whether a tensor or Lattice<Tensor> is SIMD vector or scalar
+  template<typename T, typename V=void> struct isSIMDvectorized : public std::false_type {};
+  template<typename U> struct isSIMDvectorized<U, typename std::enable_if< !std::is_same<
+    typename GridTypeMapper<typename getVectorType<U>::type>::scalar_type,
+    typename GridTypeMapper<typename getVectorType<U>::type>::vector_type>::value, void>::type>
+  : public std::true_type {};
 
-    template<typename U>
-    static double test(...);
-  
-  public:
-    enum {value = sizeof(test<T>(0)) == sizeof(char) };
-  };
-  
   //Get the precision of a Lattice, tensor or scalar type in units of sizeof(float)
   template<typename T>
   class getPrecision{

From acd25d0d01c678d1bfad72135dad8bb57a935b55 Mon Sep 17 00:00:00 2001
From: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Date: Mon, 4 Mar 2019 11:30:15 +0000
Subject: [PATCH 7/8] Wilson clover multi grid for lime lattice

---
 tests/solver/Test_wilsonclover_mg_lime.cc | 194 ++++++++++++++++++++++
 1 file changed, 194 insertions(+)
 create mode 100644 tests/solver/Test_wilsonclover_mg_lime.cc

diff --git a/tests/solver/Test_wilsonclover_mg_lime.cc b/tests/solver/Test_wilsonclover_mg_lime.cc
new file mode 100644
index 00000000..687ec83f
--- /dev/null
+++ b/tests/solver/Test_wilsonclover_mg_lime.cc
@@ -0,0 +1,194 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/solver/Test_wilsonclover_mg_mp.cc
+
+    Copyright (C) 2015-2018
+
+    Author: Daniel Richtmann <daniel.richtmann@ur.de>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+/* */
+
+
+#include <Grid/Grid.h>
+#include <Test_multigrid_common.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main(int argc, char **argv) {
+
+
+  Grid_init(&argc, &argv);
+
+  // clang-format off
+  GridCartesian         *FGrid_d   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi());
+  GridCartesian         *FGrid_f   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi());
+  GridRedBlackCartesian *FrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_d);
+  GridRedBlackCartesian *FrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid_f);
+  // clang-format on
+
+  std::vector<int> fSeeds({1, 2, 3, 4});
+  GridParallelRNG  fPRNG(FGrid_d);
+  fPRNG.SeedFixedIntegers(fSeeds);
+
+  // clang-format off
+  LatticeFermionD       src_d(FGrid_d); gaussian(fPRNG, src_d);
+  LatticeFermionD resultMGD_d(FGrid_d); resultMGD_d = zero;
+  LatticeFermionD resultMGF_d(FGrid_d); resultMGF_d = zero;
+  LatticeGaugeFieldD    Umu_d(FGrid_d);
+
+#if 0
+  {
+    FieldMetaData header;
+    std::string file("./qcdsf.769.00399.lime");
+    std::cout <<GridLogMessage<<"**************************************"<<std::endl;
+    std::cout <<GridLogMessage<<"** Reading back ILDG conf    *********"<<std::endl;
+    std::cout <<GridLogMessage<<"**************************************"<<std::endl;
+    IldgReader _IldgReader;
+    _IldgReader.open(file);
+    _IldgReader.readConfiguration(Umu_d,header);
+    _IldgReader.close();
+
+  }
+#else
+{
+  FieldMetaData header;
+  std::string file("./ckpoint_lat.IEEE64BIG.1100");
+  NerscIO::readConfiguration(Umu_d,header,file);
+}
+#endif
+  // SU3::HotConfiguration(fPRNG, Umu_d);
+
+  LatticeGaugeFieldF    Umu_f(FGrid_f); precisionChange(Umu_f, Umu_d);
+  // clang-format on
+
+  RealD mass  = -0.25;
+  RealD csw_r = 1.0;
+  RealD csw_t = 1.0;
+
+  MultiGridParams mgParams;
+  std::string     inputXml{"./mg_params.xml"};
+
+  if(GridCmdOptionExists(argv, argv + argc, "--inputxml")) {
+    inputXml = GridCmdOptionPayload(argv, argv + argc, "--inputxml");
+    assert(inputXml.length() != 0);
+  }
+
+  {
+    XmlWriter writer("mg_params_template.xml");
+    write(writer, "Params", mgParams);
+    std::cout << GridLogMessage << "Written mg_params_template.xml" << std::endl;
+
+    XmlReader reader(inputXml);
+    read(reader, "Params", mgParams);
+    std::cout << GridLogMessage << "Read in " << inputXml << std::endl;
+  }
+
+  checkParameterValidity(mgParams);
+  std::cout << mgParams << std::endl;
+
+  LevelInfo levelInfo_d(FGrid_d, mgParams);
+  LevelInfo levelInfo_f(FGrid_f, mgParams);
+
+  // Note: We do chiral doubling, so actually only nbasis/2 full basis vectors are used
+  const int nbasis = 40;
+
+  WilsonCloverFermionD Dwc_d(Umu_d, *FGrid_d, *FrbGrid_d, mass, csw_r, csw_t);
+  WilsonCloverFermionF Dwc_f(Umu_f, *FGrid_f, *FrbGrid_f, mass, csw_r, csw_t);
+
+  MdagMLinearOperator<WilsonCloverFermionD, LatticeFermionD> MdagMOpDwc_d(Dwc_d);
+  MdagMLinearOperator<WilsonCloverFermionF, LatticeFermionF> MdagMOpDwc_f(Dwc_f);
+
+  std::cout << GridLogMessage << "**************************************************" << std::endl;
+  std::cout << GridLogMessage << "Testing single-precision Multigrid for Wilson Clover" << std::endl;
+  std::cout << GridLogMessage << "**************************************************" << std::endl;
+
+  auto MGPreconDwc_f = createMGInstance<vSpinColourVectorF, vTComplexF, nbasis, WilsonCloverFermionF>(mgParams, levelInfo_f, Dwc_f, Dwc_f);
+
+  MGPreconDwc_f->setup();
+
+  if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
+    MGPreconDwc_f->runChecks(1e-6);
+  }
+
+  MixedPrecisionFlexibleGeneralisedMinimalResidual<LatticeFermionD, LatticeFermionF> MPFGMRESPREC(
+    1.0e-12, 50000, FGrid_f, *MGPreconDwc_f, 100, false);
+
+  std::cout << std::endl << "Starting with a new solver" << std::endl;
+  MPFGMRESPREC(MdagMOpDwc_d, src_d, resultMGF_d);
+
+  MGPreconDwc_f->reportTimings();
+
+  if(GridCmdOptionExists(argv, argv + argc, "--docomparison")) {
+
+    std::cout << GridLogMessage << "**************************************************" << std::endl;
+    std::cout << GridLogMessage << "Testing double-precision Multigrid for Wilson Clover" << std::endl;
+    std::cout << GridLogMessage << "**************************************************" << std::endl;
+
+    auto MGPreconDwc_d = createMGInstance<vSpinColourVectorD, vTComplexD, nbasis, WilsonCloverFermionD>(mgParams, levelInfo_d, Dwc_d, Dwc_d);
+
+    MGPreconDwc_d->setup();
+
+    if(GridCmdOptionExists(argv, argv + argc, "--runchecks")) {
+      MGPreconDwc_d->runChecks(1e-13);
+    }
+
+    FlexibleGeneralisedMinimalResidual<LatticeFermionD> FGMRESPREC(1.0e-12, 50000, *MGPreconDwc_d, 100, false);
+
+    std::cout << std::endl << "Starting with a new solver" << std::endl;
+    FGMRESPREC(MdagMOpDwc_d, src_d, resultMGD_d);
+
+    MGPreconDwc_d->reportTimings();
+
+    std::cout << GridLogMessage << "**************************************************" << std::endl;
+    std::cout << GridLogMessage << "Comparing single-precision Multigrid with double-precision one for Wilson Clover" << std::endl;
+    std::cout << GridLogMessage << "**************************************************" << std::endl;
+
+    LatticeFermionD diffFullSolver(FGrid_d);
+
+    RealD deviationFullSolver = axpy_norm(diffFullSolver, -1.0, resultMGF_d, resultMGD_d);
+
+    // clang-format off
+    LatticeFermionF src_f(FGrid_f);    precisionChange(src_f, src_d);
+    LatticeFermionF resMGF_f(FGrid_f); resMGF_f = zero;
+    LatticeFermionD resMGD_d(FGrid_d); resMGD_d = zero;
+    // clang-format on
+
+    (*MGPreconDwc_f)(src_f, resMGF_f);
+    (*MGPreconDwc_d)(src_d, resMGD_d);
+
+    LatticeFermionD diffOnlyMG(FGrid_d);
+    LatticeFermionD resMGF_d(FGrid_d);
+    precisionChange(resMGF_d, resMGF_f);
+
+    RealD deviationOnlyPrec = axpy_norm(diffOnlyMG, -1.0, resMGF_d, resMGD_d);
+
+    // clang-format off
+    std::cout << GridLogMessage << "Absolute difference between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver                      << std::endl;
+    std::cout << GridLogMessage << "Relative deviation  between FGMRES preconditioned by double and single precicision MG: " << deviationFullSolver / norm2(resultMGD_d) << std::endl;
+    std::cout << GridLogMessage << "Absolute difference between one iteration of MG Prec in double and single precision:   " << deviationOnlyPrec                        << std::endl;
+    std::cout << GridLogMessage << "Relative deviation  between one iteration of MG Prec in double and single precision:   " << deviationOnlyPrec / norm2(resMGD_d)      << std::endl;
+    // clang-format on
+  }
+
+  Grid_finalize();
+}

From 91cffef883a60f83058142404b1166e467a70d56 Mon Sep 17 00:00:00 2001
From: Michael Marshall <michael.marshall@ed.ac.uk>
Date: Thu, 7 Mar 2019 14:30:35 +0000
Subject: [PATCH 8/8] Updates after review with Peter.

---
 Grid/simd/Grid_vector_types.h |   3 +
 Grid/tensors/Tensor_class.h   | 179 ++++++----------------------------
 Grid/tensors/Tensor_traits.h  | 115 +++++++++++-----------
 3 files changed, 95 insertions(+), 202 deletions(-)

diff --git a/Grid/simd/Grid_vector_types.h b/Grid/simd/Grid_vector_types.h
index b34e056d..707af211 100644
--- a/Grid/simd/Grid_vector_types.h
+++ b/Grid/simd/Grid_vector_types.h
@@ -10,6 +10,7 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Guido Cossu <cossu@iroiro-pc.kek.jp>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: neo <cossu@post.kek.jp>
+Author: Michael Marshall <michael.marshall@ed.ac.au>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -865,8 +866,10 @@ template <typename T>
 struct is_simd : public std::false_type {};
 template <> struct is_simd<vRealF>     : public std::true_type {};
 template <> struct is_simd<vRealD>     : public std::true_type {};
+template <> struct is_simd<vRealH>     : public std::true_type {};
 template <> struct is_simd<vComplexF>  : public std::true_type {};
 template <> struct is_simd<vComplexD>  : public std::true_type {};
+template <> struct is_simd<vComplexH>  : public std::true_type {};
 template <> struct is_simd<vInteger>   : public std::true_type {};
 
 template <typename T> using IfSimd    = Invoke<std::enable_if<is_simd<T>::value, int> >;
diff --git a/Grid/tensors/Tensor_class.h b/Grid/tensors/Tensor_class.h
index c87503c6..d59640df 100644
--- a/Grid/tensors/Tensor_class.h
+++ b/Grid/tensors/Tensor_class.h
@@ -5,6 +5,7 @@ Copyright (C) 2015
 
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Michael Marshall <michael.marshall@ed.ac.au>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -42,27 +43,26 @@ namespace Grid {
 //
 class GridTensorBase {};
 
+// Too late to remove these traits from Grid Tensors, so inherit from GridTypeMapper
+#define GridVector_CopyTraits \
+  using element = vtype; \
+  using scalar_type     = typename Traits::scalar_type; \
+  using vector_type     = typename Traits::vector_type; \
+  using vector_typeD    = typename Traits::vector_typeD; \
+  using tensor_reduced  = typename Traits::tensor_reduced; \
+  using scalar_object   = typename Traits::scalar_object; \
+  using Complexified    = typename Traits::Complexified; \
+  using Realified       = typename Traits::Realified; \
+  using DoublePrecision = typename Traits::DoublePrecision; \
+  static constexpr int TensorLevel = Traits::TensorLevel
+
 template <class vtype>
 class iScalar {
  public:
   vtype _internal;
 
-  typedef vtype element;
-  typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
-  typedef typename GridTypeMapper<vtype>::vector_type vector_type;
-  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
-  typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
-  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
-  typedef iScalar<tensor_reduced_v> tensor_reduced;
-  typedef iScalar<recurse_scalar_object> scalar_object;
-  // substitutes a real or complex version with same tensor structure
-  typedef iScalar<typename GridTypeMapper<vtype>::Complexified> Complexified;
-  typedef iScalar<typename GridTypeMapper<vtype>::Realified> Realified;
-
-  // get double precision version
-  typedef iScalar<typename GridTypeMapper<vtype>::DoublePrecision> DoublePrecision;
-  
-  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
+  using Traits = GridTypeMapper<iScalar<vtype> >;
+  GridVector_CopyTraits;
 
   // Scalar no action
   //  template<int Level> using tensor_reduce_level = typename
@@ -173,37 +173,10 @@ class iScalar {
     return stream;
   };
 
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return &_internal; }
-
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return _internal.begin(); }
-
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return (&_internal) + 1; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
-
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return &_internal; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return _internal.begin(); }
-  
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return (&_internal) + 1; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return _internal.begin() + sizeof(_internal)/sizeof(scalar_type); }
+  strong_inline const scalar_type * begin() const { return reinterpret_cast<const scalar_type *>(&_internal); }
+  strong_inline       scalar_type * begin()       { return reinterpret_cast<      scalar_type *>(&_internal); }
+  strong_inline const scalar_type * end()   const { return begin() + Traits::count; }
+  strong_inline       scalar_type * end()         { return begin() + Traits::count; }
 };
 ///////////////////////////////////////////////////////////
 // Allows to turn scalar<scalar<scalar<double>>>> back to double.
@@ -224,22 +197,9 @@ class iVector {
  public:
   vtype _internal[N];
 
-  typedef vtype element;
-  typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
-  typedef typename GridTypeMapper<vtype>::vector_type vector_type;
-  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
-  typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
-  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
-  typedef iScalar<tensor_reduced_v> tensor_reduced;
-  typedef iVector<recurse_scalar_object, N> scalar_object;
+  using Traits = GridTypeMapper<iVector<vtype, N> >;
+  GridVector_CopyTraits;
 
-  // substitutes a real or complex version with same tensor structure
-  typedef iVector<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
-  typedef iVector<typename GridTypeMapper<vtype>::Realified, N> Realified;
-
-  // get double precision version
-  typedef iVector<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
-  
   template <class T, typename std::enable_if<!isGridTensor<T>::value, T>::type
                          * = nullptr>
   strong_inline auto operator=(T arg) -> iVector<vtype, N> {
@@ -248,7 +208,6 @@ class iVector {
     return *this;
   }
 
-  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
   iVector(const Zero &z) { *this = zero; };
   iVector() = default;
   /*
@@ -334,37 +293,10 @@ class iVector {
   //      return _internal[i];
   //    }
 
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return _internal; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return _internal[0].begin(); }
-  
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return _internal + N; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
-
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return _internal; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return _internal[0].begin(); }
-  
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return _internal + N; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return _internal[0].begin() + sizeof(_internal)/sizeof(scalar_type); }
+  strong_inline const scalar_type * begin() const { return reinterpret_cast<const scalar_type *>(_internal); }
+  strong_inline       scalar_type * begin()       { return reinterpret_cast<      scalar_type *>(_internal); }
+  strong_inline const scalar_type * end()   const { return begin() + Traits::count; }
+  strong_inline       scalar_type * end()         { return begin() + Traits::count; }
 };
 
 template <class vtype, int N>
@@ -372,25 +304,8 @@ class iMatrix {
  public:
   vtype _internal[N][N];
 
-  typedef vtype element;
-  typedef typename GridTypeMapper<vtype>::scalar_type scalar_type;
-  typedef typename GridTypeMapper<vtype>::vector_type vector_type;
-  typedef typename GridTypeMapper<vtype>::vector_typeD vector_typeD;
-  typedef typename GridTypeMapper<vtype>::tensor_reduced tensor_reduced_v;
-  typedef typename GridTypeMapper<vtype>::scalar_object recurse_scalar_object;
-
-  // substitutes a real or complex version with same tensor structure
-  typedef iMatrix<typename GridTypeMapper<vtype>::Complexified, N> Complexified;
-  typedef iMatrix<typename GridTypeMapper<vtype>::Realified, N> Realified;
-
-  // get double precision version
-  typedef iMatrix<typename GridTypeMapper<vtype>::DoublePrecision, N> DoublePrecision;
-  
-  // Tensor removal
-  typedef iScalar<tensor_reduced_v> tensor_reduced;
-  typedef iMatrix<recurse_scalar_object, N> scalar_object;
-
-  static constexpr int TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1;
+  using Traits = GridTypeMapper<iMatrix<vtype, N> >;
+  GridVector_CopyTraits;
 
   iMatrix(const Zero &z) { *this = zero; };
   iMatrix() = default;
@@ -521,37 +436,10 @@ class iMatrix {
   //    return _internal[i][j];
   //  }
 
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return _internal[0]; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline begin() const { return _internal[0][0].begin(); }
-  
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return _internal[0] + N * N; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, const scalar_type *>::type
-  strong_inline end() const { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
-
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return _internal[0]; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline begin() { return _internal[0][0].begin(); }
-  
-  template <typename T = vtype>
-  typename std::enable_if<!isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return _internal[0] + N * N; }
-  
-  template <typename T = vtype>
-  typename std::enable_if<isGridTensor<T>::value, scalar_type *>::type
-  strong_inline end() { return _internal[0][0].begin() + sizeof(_internal)/sizeof(scalar_type); }
+  strong_inline const scalar_type * begin() const { return reinterpret_cast<const scalar_type *>(_internal[0]); }
+  strong_inline       scalar_type * begin()       { return reinterpret_cast<      scalar_type *>(_internal[0]); }
+  strong_inline const scalar_type * end()   const { return begin() + Traits::count; }
+  strong_inline       scalar_type * end()         { return begin() + Traits::count; }
 };
 
 template <class v>
@@ -574,6 +462,3 @@ void vprefetch(const iMatrix<v, N> &vv) {
 }
 }
 #endif
-
-
-
diff --git a/Grid/tensors/Tensor_traits.h b/Grid/tensors/Tensor_traits.h
index d4854768..9cb93e17 100644
--- a/Grid/tensors/Tensor_traits.h
+++ b/Grid/tensors/Tensor_traits.h
@@ -5,6 +5,7 @@
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Christopher Kelly <ckelly@phys.columbia.edu>
+Author: Michael Marshall <michael.marshall@ed.ac.au>
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
@@ -32,14 +33,10 @@ namespace Grid {
   template<class T, int N> class iMatrix;
 
   // These are the Grid tensors
-  template<typename T> struct isGridTensor
-  : public std::false_type { static constexpr bool notvalue = true; };
-  template<class T> struct isGridTensor<iScalar<T>>
-  : public std::true_type { static constexpr bool notvalue = false; };
-  template<class T, int N> struct isGridTensor<iVector<T, N>>
-  : public std::true_type { static constexpr bool notvalue = false; };
-  template<class T, int N> struct isGridTensor<iMatrix<T, N>>
-  : public std::true_type { static constexpr bool notvalue = false; };
+  template<typename T>     struct isGridTensor                : public std::false_type { static constexpr bool notvalue = true; };
+  template<class T>        struct isGridTensor<iScalar<T>>    : public std::true_type  { static constexpr bool notvalue = false; };
+  template<class T, int N> struct isGridTensor<iVector<T, N>> : public std::true_type  { static constexpr bool notvalue = false; };
+  template<class T, int N> struct isGridTensor<iMatrix<T, N>> : public std::true_type  { static constexpr bool notvalue = false; };
 
 //////////////////////////////////////////////////////////////////////////////////
 // Want to recurse: GridTypeMapper<Matrix<vComplexD> >::scalar_type == ComplexD.
@@ -57,20 +54,15 @@ namespace Grid {
 //////////////////////////////////////////////////////////////////////////////////
 
   // This saves repeating common properties for supported Grid Scalar types
-  template<typename T, typename V=void> struct GridTypeMapper_Base {};
   // TensorLevel    How many nested grid tensors
   // Rank           Rank of the grid tensor
   // count          Total number of elements, i.e. product of dimensions
-  // scalar_size    Size of the underlying fundamental object (tensor_reduced) in bytes
-  // size           Total size of all elements in bytes
   // Dimension(dim) Size of dimension dim
-  template<typename T> struct GridTypeMapper_Base<T> {
+  struct GridTypeMapper_Base {
     static constexpr int TensorLevel = 0;
     static constexpr int Rank = 0;
     static constexpr std::size_t count = 1;
-    static constexpr std::size_t scalar_size = sizeof(T);
-    static constexpr std::size_t size = scalar_size * count;
-    static constexpr int Dimension(unsigned int dim) { return 0; }
+    static constexpr int Dimension(int dim) { return 0; }
   };
 
 //////////////////////////////////////////////////////////////////////////////////
@@ -79,7 +71,7 @@ namespace Grid {
 
   template<typename T> struct GridTypeMapper {};
 
-  template<> struct GridTypeMapper<RealF> : public GridTypeMapper_Base<RealF> {
+  template<> struct GridTypeMapper<RealF> : public GridTypeMapper_Base {
     typedef RealF scalar_type;
     typedef RealF vector_type;
     typedef RealD vector_typeD;
@@ -89,7 +81,7 @@ namespace Grid {
     typedef RealF Realified;
     typedef RealD DoublePrecision;
   };
-  template<> struct GridTypeMapper<RealD> : public GridTypeMapper_Base<RealD> {
+  template<> struct GridTypeMapper<RealD> : public GridTypeMapper_Base {
     typedef RealD scalar_type;
     typedef RealD vector_type;
     typedef RealD vector_typeD;
@@ -99,7 +91,7 @@ namespace Grid {
     typedef RealD Realified;
     typedef RealD DoublePrecision;
   };
-  template<> struct GridTypeMapper<ComplexF> : public GridTypeMapper_Base<ComplexF> {
+  template<> struct GridTypeMapper<ComplexF> : public GridTypeMapper_Base {
     typedef ComplexF scalar_type;
     typedef ComplexF vector_type;
     typedef ComplexD vector_typeD;
@@ -109,7 +101,7 @@ namespace Grid {
     typedef RealF Realified;
     typedef ComplexD DoublePrecision;
   };
-  template<> struct GridTypeMapper<ComplexD> : public GridTypeMapper_Base<ComplexD> {
+  template<> struct GridTypeMapper<ComplexD> : public GridTypeMapper_Base {
     typedef ComplexD scalar_type;
     typedef ComplexD vector_type;
     typedef ComplexD vector_typeD;
@@ -119,7 +111,7 @@ namespace Grid {
     typedef RealD Realified;
     typedef ComplexD DoublePrecision;
   };
-  template<> struct GridTypeMapper<Integer> : public GridTypeMapper_Base<Integer> {
+  template<> struct GridTypeMapper<Integer> : public GridTypeMapper_Base {
     typedef Integer scalar_type;
     typedef Integer vector_type;
     typedef Integer vector_typeD;
@@ -130,7 +122,7 @@ namespace Grid {
     typedef void DoublePrecision;
   };
 
-  template<> struct GridTypeMapper<vRealF> : public GridTypeMapper_Base<vRealF> {
+  template<> struct GridTypeMapper<vRealF> : public GridTypeMapper_Base {
     typedef RealF  scalar_type;
     typedef vRealF vector_type;
     typedef vRealD vector_typeD;
@@ -140,7 +132,7 @@ namespace Grid {
     typedef vRealF Realified;
     typedef vRealD DoublePrecision;
   };
-  template<> struct GridTypeMapper<vRealD> : public GridTypeMapper_Base<vRealD> {
+  template<> struct GridTypeMapper<vRealD> : public GridTypeMapper_Base {
     typedef RealD  scalar_type;
     typedef vRealD vector_type;
     typedef vRealD vector_typeD;
@@ -150,7 +142,17 @@ namespace Grid {
     typedef vRealD Realified;
     typedef vRealD DoublePrecision;
   };
-  template<> struct GridTypeMapper<vComplexH> : public GridTypeMapper_Base<vComplexH> {
+  template<> struct GridTypeMapper<vRealH> : public GridTypeMapper_Base {
+    typedef RealF  scalar_type;
+    typedef vRealH vector_type;
+    typedef vRealD vector_typeD;
+    typedef vRealH tensor_reduced;
+    typedef RealF  scalar_object;
+    typedef vComplexH Complexified;
+    typedef vRealH Realified;
+    typedef vRealD DoublePrecision;
+  };
+  template<> struct GridTypeMapper<vComplexH> : public GridTypeMapper_Base {
     typedef ComplexF  scalar_type;
     typedef vComplexH vector_type;
     typedef vComplexD vector_typeD;
@@ -160,7 +162,7 @@ namespace Grid {
     typedef vRealH Realified;
     typedef vComplexD DoublePrecision;
   };
-  template<> struct GridTypeMapper<vComplexF> : public GridTypeMapper_Base<vComplexF> {
+  template<> struct GridTypeMapper<vComplexF> : public GridTypeMapper_Base {
     typedef ComplexF  scalar_type;
     typedef vComplexF vector_type;
     typedef vComplexD vector_typeD;
@@ -170,7 +172,7 @@ namespace Grid {
     typedef vRealF Realified;
     typedef vComplexD DoublePrecision;
   };
-  template<> struct GridTypeMapper<vComplexD> : public GridTypeMapper_Base<vComplexD> {
+  template<> struct GridTypeMapper<vComplexD> : public GridTypeMapper_Base {
     typedef ComplexD  scalar_type;
     typedef vComplexD vector_type;
     typedef vComplexD vector_typeD;
@@ -180,7 +182,7 @@ namespace Grid {
     typedef vRealD Realified;
     typedef vComplexD DoublePrecision;
   };
-  template<> struct GridTypeMapper<vInteger> : public GridTypeMapper_Base<vInteger> {
+  template<> struct GridTypeMapper<vInteger> : public GridTypeMapper_Base {
     typedef  Integer scalar_type;
     typedef vInteger vector_type;
     typedef vInteger vector_typeD;
@@ -192,46 +194,49 @@ namespace Grid {
   };
 
 #define GridTypeMapper_RepeatedTypes \
-  typedef typename ObjectTraits::scalar_type scalar_type; \
-  typedef typename ObjectTraits::vector_type vector_type; \
-  typedef typename ObjectTraits::vector_typeD vector_typeD; \
-  typedef typename ObjectTraits::tensor_reduced tensor_reduced; \
-  typedef typename ObjectTraits::scalar_object scalar_object; \
-  typedef typename ObjectTraits::Complexified Complexified; \
-  typedef typename ObjectTraits::Realified Realified; \
-  typedef typename ObjectTraits::DoublePrecision DoublePrecision; \
-  static constexpr int TensorLevel = BaseTraits::TensorLevel + 1; \
-  static constexpr std::size_t scalar_size = BaseTraits::scalar_size; \
-  static constexpr std::size_t size = scalar_size * count
+  using BaseTraits   = GridTypeMapper<T>; \
+  using scalar_type  = typename BaseTraits::scalar_type; \
+  using vector_type  = typename BaseTraits::vector_type; \
+  using vector_typeD = typename BaseTraits::vector_typeD; \
+  static constexpr int TensorLevel = BaseTraits::TensorLevel + 1
 
   template<typename T> struct GridTypeMapper<iScalar<T>> {
-    using ObjectTraits = iScalar<T>;
-    using BaseTraits = GridTypeMapper<T>;
-    static constexpr int Rank = 1 + BaseTraits::Rank;
-    static constexpr std::size_t count = 1 * BaseTraits::count;
-    static constexpr int Dimension(unsigned int dim) {
-      return ( dim == 0 ) ? 1 : BaseTraits::Dimension(dim - 1); }
     GridTypeMapper_RepeatedTypes;
+    using tensor_reduced  = iScalar<typename BaseTraits::tensor_reduced>;
+    using scalar_object   = iScalar<typename BaseTraits::scalar_object>;
+    using Complexified    = iScalar<typename BaseTraits::Complexified>;
+    using Realified       = iScalar<typename BaseTraits::Realified>;
+    using DoublePrecision = iScalar<typename BaseTraits::DoublePrecision>;
+    static constexpr int Rank = BaseTraits::Rank + 1;
+    static constexpr std::size_t count = BaseTraits::count;
+    static constexpr int Dimension(int dim) {
+      return ( dim == 0 ) ? 1 : BaseTraits::Dimension(dim - 1); }
   };
 
   template<typename T, int N> struct GridTypeMapper<iVector<T, N>> {
-    using ObjectTraits = iVector<T, N>;
-    using BaseTraits = GridTypeMapper<T>;
-    static constexpr int Rank = 1 + BaseTraits::Rank;
-    static constexpr std::size_t count = N * BaseTraits::count;
-    static constexpr int Dimension(unsigned int dim) {
-      return ( dim == 0 ) ? N : BaseTraits::Dimension(dim - 1); }
     GridTypeMapper_RepeatedTypes;
+    using tensor_reduced  = iScalar<typename BaseTraits::tensor_reduced>;
+    using scalar_object   = iVector<typename BaseTraits::scalar_object,   N>;
+    using Complexified    = iVector<typename BaseTraits::Complexified,    N>;
+    using Realified       = iVector<typename BaseTraits::Realified,       N>;
+    using DoublePrecision = iVector<typename BaseTraits::DoublePrecision, N>;
+    static constexpr int Rank = BaseTraits::Rank + 1;
+    static constexpr std::size_t count = BaseTraits::count * N;
+    static constexpr int Dimension(int dim) {
+      return ( dim == 0 ) ? N : BaseTraits::Dimension(dim - 1); }
   };
 
   template<typename T, int N> struct GridTypeMapper<iMatrix<T, N>> {
-    using ObjectTraits = iMatrix<T, N>;
-    using BaseTraits = GridTypeMapper<T>;
-    static constexpr int Rank = 2 + BaseTraits::Rank;
-    static constexpr std::size_t count = N * N * BaseTraits::count;
-    static constexpr int Dimension(unsigned int dim) {
-      return ( dim == 0 || dim == 1 ) ? N : BaseTraits::Dimension(dim - 2); }
     GridTypeMapper_RepeatedTypes;
+    using tensor_reduced  = iScalar<typename BaseTraits::tensor_reduced>;
+    using scalar_object   = iMatrix<typename BaseTraits::scalar_object,   N>;
+    using Complexified    = iMatrix<typename BaseTraits::Complexified,    N>;
+    using Realified       = iMatrix<typename BaseTraits::Realified,       N>;
+    using DoublePrecision = iMatrix<typename BaseTraits::DoublePrecision, N>;
+    static constexpr int Rank = BaseTraits::Rank + 2;
+    static constexpr std::size_t count = BaseTraits::count * N * N;
+    static constexpr int Dimension(int dim) {
+      return ( dim == 0 || dim == 1 ) ? N : BaseTraits::Dimension(dim - 2); }
   };
 
   // Match the index