Merge branch 'master' of https://github.com/paboyle/Grid

2025-11-04 14:04:32 +00:00 · 2015-06-05 10:04:46 +01:00
parent 7d984b9547 f22170ad49
commit db84b19443
150 changed files with 5950 additions and 2966 deletions
--- a/lib/.dirstamp
+++ b/lib/.dirstamp
--- a/lib/Grid_algorithms.h
+++ b/lib/Grid_algorithms.h
--- a/lib/Grid_aligned_allocator.h
+++ b/lib/Grid_aligned_allocator.h
--- a/lib/Cartesian.h
+++ b/lib/Cartesian.h
@@ -0,0 +1,8 @@
+#ifndef GRID_CARTESIAN_H
+#define GRID_CARTESIAN_H
+
+#include <cartesian/Cartesian_base.h>
+#include <cartesian/Cartesian_full.h>
+#include <cartesian/Cartesian_red_black.h> 
+
+#endif
--- a/lib/Grid_communicator.h
+++ b/lib/Grid_communicator.h
@@ -1,6 +1,6 @@
 #ifndef GRID_COMMUNICATOR_H
 #define GRID_COMMUNICATOR_H

-#include <communicator/Grid_communicator_base.h>
+#include <communicator/Communicator_base.h>

 #endif
--- a/lib/Grid_comparison.h
+++ b/lib/Grid_comparison.h
@@ -141,7 +141,7 @@ namespace Grid {
    }
 }

-#include <lattice/Grid_lattice_comparison.h>
-#include <lattice/Grid_lattice_where.h>
+#include <lattice/Lattice_comparison.h>
+#include <lattice/Lattice_where.h>

 #endif
--- a/lib/Grid_cshift.h
+++ b/lib/Grid_cshift.h
@@ -1,13 +1,13 @@
 #ifndef _GRID_CSHIFT_H_
 #define _GRID_CSHIFT_H_

-#include <cshift/Grid_cshift_common.h>
+#include <cshift/Cshift_common.h>

 #ifdef GRID_COMMS_NONE
-#include <cshift/Grid_cshift_none.h>
+#include <cshift/Cshift_none.h>
 #endif

 #ifdef GRID_COMMS_MPI
-#include <cshift/Grid_cshift_mpi.h>
+#include <cshift/Cshift_mpi.h>
 #endif 
 #endif
--- a/lib/Grid.h
+++ b/lib/Grid.h
@@ -33,7 +33,7 @@

 #define strong_inline __attribute__((always_inline)) inline

-#include <Grid_config.h>
+#include <GridConfig.h>

 ////////////////////////////////////////////////////////////
 // Tunable header includes
@@ -46,22 +46,22 @@
 #include <malloc.h>
 #endif

-#include <Grid_aligned_allocator.h>
+#include <AlignedAllocator.h>

-#include <Grid_simd.h>
-#include <Grid_threads.h>
+#include <Simd.h>
+#include <Threads.h>

-#include <Grid_cartesian.h> // subdir aggregate
-#include <Grid_math.h>      // subdir aggregate
-#include <Grid_lattice.h>   // subdir aggregate
-#include <Grid_comparison.h>
-#include <Grid_cshift.h>    // subdir aggregate
-#include <Grid_stencil.h>   // subdir aggregate
-
-#include <Grid_algorithms.h>// subdir aggregate
+#include <Communicator.h> // subdir aggregate
+#include <Cartesian.h> // subdir aggregate
+#include <Tensors.h>   // subdir aggregate
+#include <Lattice.h>   // subdir aggregate
+#include <Comparison.h>
+#include <Cshift.h>    // subdir aggregate
+#include <Stencil.h>   // subdir aggregate
+#include <Algorithms.h>// subdir aggregate

 #include <qcd/QCD.h>
-#include <parallelIO/GridNerscIO.h>
+#include <parallelIO/NerscIO.h>

 namespace Grid {

--- a/lib/Grid_config.h
+++ b/lib/Grid_config.h
@@ -1,5 +1,5 @@
-/* lib/Grid_config.h.  Generated from Grid_config.h.in by configure.  */
-/* lib/Grid_config.h.in.  Generated from configure.ac by autoheader.  */
+/* lib/GridConfig.h.  Generated from GridConfig.h.in by configure.  */
+/* lib/GridConfig.h.in.  Generated from configure.ac by autoheader.  */

 /* AVX */
 /* #undef AVX1 */
@@ -16,6 +16,15 @@
 /* GRID_COMMS_NONE */
 #define GRID_COMMS_NONE 1

+/* Support Altivec instructions */
+/* #undef HAVE_ALTIVEC */
+
+/* Support AVX (Advanced Vector Extensions) instructions */
+/* #undef HAVE_AVX */
+
+/* Support AVX2 (Advanced Vector Extensions 2) instructions */
+/* #undef HAVE_AVX2 */
+
 /* define if the compiler supports basic C++11 syntax */
 /* #undef HAVE_CXX11 */

@@ -30,6 +39,9 @@
 /* Define to 1 if you have the <endian.h> header file. */
 #define HAVE_ENDIAN_H 1

+/* Support FMA3 (Fused Multiply-Add) instructions */
+/* #undef HAVE_FMA */
+
 /* Define to 1 if you have the `gettimeofday' function. */
 #define HAVE_GETTIMEOFDAY 1

@@ -54,9 +66,30 @@
 /* Define to 1 if you have the <memory.h> header file. */
 #define HAVE_MEMORY_H 1

+/* Support mmx instructions */
+#define HAVE_MMX /**/
+
 /* Define to 1 if you have the <mm_malloc.h> header file. */
 #define HAVE_MM_MALLOC_H 1

+/* Support SSE (Streaming SIMD Extensions) instructions */
+#define HAVE_SSE /**/
+
+/* Support SSE2 (Streaming SIMD Extensions 2) instructions */
+#define HAVE_SSE2 /**/
+
+/* Support SSE3 (Streaming SIMD Extensions 3) instructions */
+#define HAVE_SSE3 /**/
+
+/* Support SSSE4.1 (Streaming SIMD Extensions 4.1) instructions */
+#define HAVE_SSE4_1 /**/
+
+/* Support SSSE4.2 (Streaming SIMD Extensions 4.2) instructions */
+#define HAVE_SSE4_2 /**/
+
+/* Support SSSE3 (Supplemental Streaming SIMD Extensions 3) instructions */
+#define HAVE_SSSE3 /**/
+
 /* Define to 1 if you have the <stdint.h> header file. */
 #define HAVE_STDINT_H 1

--- a/lib/Grid_config.h.in
+++ b/lib/Grid_config.h.in
@@ -1,4 +1,4 @@
-/* lib/Grid_config.h.in.  Generated from configure.ac by autoheader.  */
+/* lib/GridConfig.h.in.  Generated from configure.ac by autoheader.  */

 /* AVX */
 #undef AVX1
@@ -15,6 +15,15 @@
 /* GRID_COMMS_NONE */
 #undef GRID_COMMS_NONE

+/* Support Altivec instructions */
+#undef HAVE_ALTIVEC
+
+/* Support AVX (Advanced Vector Extensions) instructions */
+#undef HAVE_AVX
+
+/* Support AVX2 (Advanced Vector Extensions 2) instructions */
+#undef HAVE_AVX2
+
 /* define if the compiler supports basic C++11 syntax */
 #undef HAVE_CXX11

@@ -29,6 +38,9 @@
 /* Define to 1 if you have the <endian.h> header file. */
 #undef HAVE_ENDIAN_H

+/* Support FMA3 (Fused Multiply-Add) instructions */
+#undef HAVE_FMA
+
 /* Define to 1 if you have the `gettimeofday' function. */
 #undef HAVE_GETTIMEOFDAY

@@ -53,9 +65,30 @@
 /* Define to 1 if you have the <memory.h> header file. */
 #undef HAVE_MEMORY_H

+/* Support mmx instructions */
+#undef HAVE_MMX
+
 /* Define to 1 if you have the <mm_malloc.h> header file. */
 #undef HAVE_MM_MALLOC_H

+/* Support SSE (Streaming SIMD Extensions) instructions */
+#undef HAVE_SSE
+
+/* Support SSE2 (Streaming SIMD Extensions 2) instructions */
+#undef HAVE_SSE2
+
+/* Support SSE3 (Streaming SIMD Extensions 3) instructions */
+#undef HAVE_SSE3
+
+/* Support SSSE4.1 (Streaming SIMD Extensions 4.1) instructions */
+#undef HAVE_SSE4_1
+
+/* Support SSSE4.2 (Streaming SIMD Extensions 4.2) instructions */
+#undef HAVE_SSE4_2
+
+/* Support SSSE3 (Supplemental Streaming SIMD Extensions 3) instructions */
+#undef HAVE_SSSE3
+
 /* Define to 1 if you have the <stdint.h> header file. */
 #undef HAVE_STDINT_H

--- a/lib/Grid_init.cc
+++ b/lib/Grid_init.cc
@@ -143,7 +143,7 @@ void Grid_init(int *argc,char ***argv)
  }
  if( GridCmdOptionExists(*argv,*argv+*argc,"--dslash-opt") ){
    WilsonFermion::HandOptDslash=1;
-    FiveDimWilsonFermion::HandOptDslash=1;
+    WilsonFermion5D::HandOptDslash=1;
  }
  if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){
    LebesgueOrder::UseLebesgueOrder=1;
--- a/lib/Grid_cartesian.h
+++ b/lib/Grid_cartesian.h
@@ -1,8 +0,0 @@
-#ifndef GRID_CARTESIAN_H
-#define GRID_CARTESIAN_H
-
-#include <cartesian/Grid_cartesian_base.h>
-#include <cartesian/Grid_cartesian_full.h>
-#include <cartesian/Grid_cartesian_red_black.h> 
-
-#endif
--- a/lib/Grid_math.h
+++ b/lib/Grid_math.h
@@ -1,16 +0,0 @@
-#ifndef GRID_MATH_H
-#define GRID_MATH_H
-
-#include <math/Grid_math_traits.h>
-#include <math/Grid_math_tensors.h>
-#include <math/Grid_math_arith.h>
-#include <math/Grid_math_inner.h>
-#include <math/Grid_math_outer.h>
-#include <math/Grid_math_transpose.h>
-#include <math/Grid_math_trace.h>
-#include <math/Grid_math_peek.h>
-#include <math/Grid_math_poke.h>
-#include <math/Grid_math_reality.h>
-
-    
-#endif
--- a/lib/Grid_lattice.h
+++ b/lib/Grid_lattice.h
@@ -1,6 +1,6 @@
 #ifndef GRID_LATTICE_H
 #define GRID_LATTICE_H

-#include <lattice/Grid_lattice_base.h>
+#include <lattice/Lattice_base.h>

 #endif
--- a/lib/Make.inc
+++ b/lib/Make.inc
@@ -0,0 +1,4 @@
+
+HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./Comparison.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./GridConfig.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_where.h ./Lattice.h ./parallelIO/NerscIO.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/Dirac.h ./qcd/LinalgUtils.h ./qcd/QCD.h ./qcd/SpaceTimeGrid.h ./qcd/TwoSpinor.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Old/Grid_vComplexD.h ./simd/Old/Grid_vComplexF.h ./simd/Old/Grid_vInteger.h ./simd/Old/Grid_vRealD.h ./simd/Old/Grid_vRealF.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_inner.h ./tensors/Tensor_outer.h ./tensors/Tensor_peek.h ./tensors/Tensor_poke.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./Tensors.h ./Threads.h
+
+CCFILES=./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/Dirac.cc ./qcd/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -3,115 +3,26 @@ AM_CXXFLAGS = -I$(top_srcdir)/

 extra_sources=
 if BUILD_COMMS_MPI
-  extra_sources+=communicator/Grid_communicator_mpi.cc
+  extra_sources+=communicator/Communicator_mpi.cc
 endif

 if BUILD_COMMS_NONE
-  extra_sources+=communicator/Grid_communicator_none.cc
+  extra_sources+=communicator/Communicator_none.cc
 endif

 #
 # Libraries
 #
-lib_LIBRARIES = libGrid.a
-libGrid_a_SOURCES =				\
-	Grid_init.cc				\
-	stencil/Grid_lebesgue.cc		\
-	stencil/Grid_stencil_common.cc		\
-	algorithms/approx/Zolotarev.cc		\
-	algorithms/approx/Remez.cc		\
-	qcd/action/fermion/FiveDimWilsonFermion.cc\
-	qcd/action/fermion/WilsonFermion.cc\
-	qcd/action/fermion/WilsonKernels.cc\
-	qcd/action/fermion/WilsonKernelsHand.cc\
-	qcd/Dirac.cc\
-	$(extra_sources)

+include Make.inc
+
+lib_LIBRARIES = libGrid.a
+libGrid_a_SOURCES = $(CCFILES) $(extra_sources)
+
+
+#	qcd/action/fermion/PartialFractionFermion5D.cc\	\
 #
 # Include files
 #
-nobase_include_HEADERS=\
-		./algorithms/approx/bigfloat.h\
-		./algorithms/approx/bigfloat_double.h\
-		./algorithms/approx/Chebyshev.h\
-		./algorithms/approx/Remez.h\
-		./algorithms/approx/Zolotarev.h\
-		./algorithms/iterative/ConjugateGradient.h\
-		./algorithms/iterative/NormalEquations.h\
-		./algorithms/iterative/SchurRedBlack.h\
-		./algorithms/LinearOperator.h\
-		./algorithms/SparseMatrix.h\
-		./cartesian/Grid_cartesian_base.h\
-		./cartesian/Grid_cartesian_full.h\
-		./cartesian/Grid_cartesian_red_black.h\
-		./communicator/Grid_communicator_base.h\
-		./cshift/Grid_cshift_common.h\
-		./cshift/Grid_cshift_mpi.h\
-		./cshift/Grid_cshift_none.h\
-		./Grid.h\
-		./Grid_algorithms.h\
-		./Grid_aligned_allocator.h\
-		./Grid_cartesian.h\
-		./Grid_communicator.h\
-		./Grid_comparison.h\
-		./Grid_config.h\
-		./Grid_cshift.h\
-		./Grid_extract.h\
-		./Grid_lattice.h\
-		./Grid_math.h\
-		./Grid_simd.h\
-		./Grid_stencil.h\
-		./Grid_threads.h\
-		./lattice/Grid_lattice_arith.h\
-		./lattice/Grid_lattice_base.h\
-		./lattice/Grid_lattice_comparison.h\
-		./lattice/Grid_lattice_conformable.h\
-		./lattice/Grid_lattice_coordinate.h\
-		./lattice/Grid_lattice_ET.h\
-		./lattice/Grid_lattice_local.h\
-		./lattice/Grid_lattice_overload.h\
-		./lattice/Grid_lattice_peekpoke.h\
-		./lattice/Grid_lattice_reality.h\
-		./lattice/Grid_lattice_reduction.h\
-		./lattice/Grid_lattice_rng.h\
-		./lattice/Grid_lattice_trace.h\
-		./lattice/Grid_lattice_transfer.h\
-		./lattice/Grid_lattice_transpose.h\
-		./lattice/Grid_lattice_where.h\
-		./math/Grid_math_arith.h\
-		./math/Grid_math_arith_add.h\
-		./math/Grid_math_arith_mac.h\
-		./math/Grid_math_arith_mul.h\
-		./math/Grid_math_arith_scalar.h\
-		./math/Grid_math_arith_sub.h\
-		./math/Grid_math_inner.h\
-		./math/Grid_math_outer.h\
-		./math/Grid_math_peek.h\
-		./math/Grid_math_poke.h\
-		./math/Grid_math_reality.h\
-		./math/Grid_math_tensors.h\
-		./math/Grid_math_trace.h\
-		./math/Grid_math_traits.h\
-		./math/Grid_math_transpose.h\
-		./parallelIO/GridNerscIO.h\
-		./qcd/action/Actions.h\
-		./qcd/action/fermion/FermionAction.h\
-		./qcd/action/fermion/FiveDimWilsonFermion.h\
-		./qcd/action/fermion/WilsonCompressor.h\
-		./qcd/action/fermion/WilsonFermion.h\
-		./qcd/action/fermion/WilsonKernels.h\
-		./qcd/Dirac.h\
-		./qcd/QCD.h\
-		./qcd/TwoSpinor.h\
-		./simd/Grid_avx.h\
-		./simd/Grid_avx512.h\
-		./simd/Grid_qpx.h\
-		./simd/Grid_sse4.h\
-		./simd/Grid_vector_types.h\
-		./simd/Old/Grid_vComplexD.h\
-		./simd/Old/Grid_vComplexF.h\
-		./simd/Old/Grid_vInteger.h\
-		./simd/Old/Grid_vRealD.h\
-		./simd/Old/Grid_vRealF.h\
-		./stencil/Grid_lebesgue.h
+nobase_include_HEADERS=$(HFILES)

--- a/lib/Grid_simd.h
+++ b/lib/Grid_simd.h
--- a/lib/Grid_stencil.h
+++ b/lib/Grid_stencil.h
@@ -1,7 +1,7 @@
 #ifndef GRID_STENCIL_H
 #define GRID_STENCIL_H

-#include <stencil/Grid_lebesgue.h>   // subdir aggregate
+#include <stencil/Lebesgue.h>   // subdir aggregate

 //////////////////////////////////////////////////////////////////////////////////////////
 // Must not lose sight that goal is to be able to construct really efficient
--- a/lib/Tensors.h
+++ b/lib/Tensors.h
@@ -0,0 +1,17 @@
+#ifndef GRID_MATH_H
+#define GRID_MATH_H
+
+#include <tensors/Tensor_traits.h>
+#include <tensors/Tensor_class.h>
+#include <tensors/Tensor_arith.h>
+#include <tensors/Tensor_inner.h>
+#include <tensors/Tensor_outer.h>
+#include <tensors/Tensor_transpose.h>
+#include <tensors/Tensor_trace.h>
+#include <tensors/Tensor_Ta.h>
+#include <tensors/Tensor_peek.h>
+#include <tensors/Tensor_poke.h>
+#include <tensors/Tensor_reality.h>
+#include <tensors/Tensor_extract_merge.h>
+    
+#endif
--- a/lib/Grid_threads.h
+++ b/lib/Grid_threads.h
--- a/lib/algorithms/LinearOperator.h
+++ b/lib/algorithms/LinearOperator.h
@@ -125,39 +125,7 @@ namespace Grid {
     };
    */

-    // Chroma interface defining GaugeAction
-    /*
-      template<typename P, typename Q>   class GaugeAction
-  virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
-  virtual GaugeState<P,Q>* createState(const Q& q) const
-  virtual const GaugeBC<P,Q>& getGaugeBC() const
-  virtual const Set& getSet(void) const = 0;
-  virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const 
-  virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;

-  class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
-  typedef multi1d<LatticeColorMatrix>  P;
-  typedef multi1d<LatticeColorMatrix>  Q;
-  virtual void staple(LatticeColorMatrix& result,
-		      const Handle< GaugeState<P,Q> >& state,
-		      int mu, int cb) const = 0;
-    */
-
-    // Chroma interface defining FermionAction
-    /*
-     template<typename T, typename P, typename Q>  class FermAct4D : public FermionAction<T,P,Q>
-     virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
-     virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
-     virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
-     virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
-     class DiffFermAct4D : public FermAct4D<T,P,Q>
-     virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
-     virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
-    */
 }

 #endif
--- a/lib/algorithms/approx/Zolotarev.cc
+++ b/lib/algorithms/approx/Zolotarev.cc
@@ -58,6 +58,8 @@

 /* Compute the partial fraction expansion coefficients (alpha) from the
 * factored form */
+namespace Grid {
+namespace Approx {

 static void construct_partfrac(izd *z) {
  int dn = z -> dn, dd = z -> dd, type = z -> type;
@@ -291,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k,
 * Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and
 * type = 1 for the approximation which is infinite at x = 0. */

-zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
+zolotarev_data* zolotarev(PRECISION epsilon, int n, int type) {
  INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F,
    l, invlambda, xi, xisq, *tv, s, opl;
  int m, czero, ts;
@@ -412,7 +414,19 @@ zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
  return zd;
 }

-zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
+
+void zolotarev_free(zolotarev_data *zdata)
+{
+    free(zdata -> a);
+    free(zdata -> ap);
+    free(zdata -> alpha);
+    free(zdata -> beta);
+    free(zdata -> gamma);
+    free(zdata);
+}
+
+
+zolotarev_data* higham(PRECISION epsilon, int n) {
  INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq;
  int m, czero;
  zolotarev_data *zd;
@@ -502,6 +516,7 @@ zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
  free(d);
  return zd;
 }
+}}

 #ifdef TEST

@@ -707,4 +722,6 @@ int main(int argc, char** argv) {

  return EXIT_SUCCESS;
 }
+
+
 #endif /* TEST */
--- a/lib/algorithms/approx/Zolotarev.h
+++ b/lib/algorithms/approx/Zolotarev.h
@@ -1,7 +1,8 @@
 /* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */

 #ifdef __cplusplus
-extern "C" {
+namespace Grid {
+namespace Approx {
 #endif

 #define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
@@ -76,10 +77,11 @@ typedef struct {
 * zolotarev_data structure. The arguments must satisfy the constraints that
 * epsilon > 0, n > 0, and type = 0 or 1. */

-ZOLOTAREV_DATA* bfm_higham(PRECISION epsilon, int n) ;
-ZOLOTAREV_DATA* bfm_zolotarev(PRECISION epsilon, int n, int type);
+ZOLOTAREV_DATA* higham(PRECISION epsilon, int n) ;
+ZOLOTAREV_DATA* zolotarev(PRECISION epsilon, int n, int type);
+void zolotarev_free(zolotarev_data *zdata);
 #endif

 #ifdef __cplusplus
-}
+}}
 #endif
--- a/lib/algorithms/iterative/NormalEquations.h
+++ b/lib/algorithms/iterative/NormalEquations.h
@@ -4,8 +4,7 @@
 namespace Grid {

  ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  // Take a matrix and form a Red Black solver calling a Herm solver
-  // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
+  // Take a matrix and form an NE solver calling a Herm solver
  ///////////////////////////////////////////////////////////////////////////////////////////////////////
  template<class Field> class NormalEquations : public OperatorFunction<Field>{
  private:
@@ -15,7 +14,7 @@ namespace Grid {
  public:

    /////////////////////////////////////////////////////
-    // Wrap the usual normal equations Schur trick
+    // Wrap the usual normal equations trick
    /////////////////////////////////////////////////////
  NormalEquations(SparseMatrixBase<Field> &Matrix, OperatorFunction<Field> &HermitianSolver) 
    :  _Matrix(Matrix), _HermitianSolver(HermitianSolver) {}; 
--- a/lib/cartesian/Grid_cartesian_base.h
+++ b/lib/cartesian/Grid_cartesian_base.h
@@ -2,7 +2,6 @@
 #define GRID_CARTESIAN_BASE_H

 #include <Grid.h>
-#include <Grid_communicator.h>

 namespace Grid{

@@ -21,7 +20,7 @@ public:
    // Give Lattice access
    template<class object> friend class Lattice;

-    GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
+    GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};


    // Physics Grid information.
--- a/lib/cartesian/Grid_cartesian_full.h
+++ b/lib/cartesian/Grid_cartesian_full.h
@@ -27,9 +27,9 @@ public:
    virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){
      return shift;
    }
-    GridCartesian(std::vector<int> &dimensions,
-		  std::vector<int> &simd_layout,
-		  std::vector<int> &processor_grid
+    GridCartesian(const std::vector<int> &dimensions,
+		  const std::vector<int> &simd_layout,
+		  const std::vector<int> &processor_grid
 		  ) : GridBase(processor_grid)
    {
        ///////////////////////
--- a/lib/cartesian/Grid_cartesian_red_black.h
+++ b/lib/cartesian/Grid_cartesian_red_black.h
@@ -81,28 +81,28 @@ public:
      }
    };

-    GridRedBlackCartesian(GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors)  {};
+    GridRedBlackCartesian(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors)  {};

-    GridRedBlackCartesian(std::vector<int> &dimensions,
-			  std::vector<int> &simd_layout,
-			  std::vector<int> &processor_grid,
-			  std::vector<int> &checker_dim_mask,
+    GridRedBlackCartesian(const std::vector<int> &dimensions,
+			  const std::vector<int> &simd_layout,
+			  const std::vector<int> &processor_grid,
+			  const std::vector<int> &checker_dim_mask,
 			  int checker_dim
 			  ) :  GridBase(processor_grid) 
    {
      Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim);
    }
-    GridRedBlackCartesian(std::vector<int> &dimensions,
-			  std::vector<int> &simd_layout,
-			  std::vector<int> &processor_grid) : GridBase(processor_grid) 
+    GridRedBlackCartesian(const std::vector<int> &dimensions,
+			  const std::vector<int> &simd_layout,
+			  const std::vector<int> &processor_grid) : GridBase(processor_grid) 
    {
      std::vector<int> checker_dim_mask(dimensions.size(),1);
      Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0);
    }
-    void Init(std::vector<int> &dimensions,
-	      std::vector<int> &simd_layout,
-	      std::vector<int> &processor_grid,
-	      std::vector<int> &checker_dim_mask,
+    void Init(const std::vector<int> &dimensions,
+	      const std::vector<int> &simd_layout,
+	      const std::vector<int> &processor_grid,
+	      const std::vector<int> &checker_dim_mask,
 	      int checker_dim)
    {
    ///////////////////////
--- a/lib/communicator/Grid_communicator_base.h
+++ b/lib/communicator/Grid_communicator_base.h
@@ -27,7 +27,7 @@ class CartesianCommunicator {
 #endif

    // Constructor
-    CartesianCommunicator(std::vector<int> &pdimensions_in);
+    CartesianCommunicator(const std::vector<int> &pdimensions_in);

    // Wraps MPI_Cart routines
    void ShiftedRanks(int dim,int shift,int & source, int & dest);
--- a/lib/communicator/Grid_communicator_mpi.cc
+++ b/lib/communicator/Grid_communicator_mpi.cc
@@ -5,7 +5,7 @@ namespace Grid {

  // Should error check all MPI calls.

-CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
+CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
 {
  _ndimension = processors.size();
  std::vector<int> periodic(_ndimension,1);
--- a/lib/communicator/Grid_communicator_none.cc
+++ b/lib/communicator/Grid_communicator_none.cc
@@ -1,7 +1,7 @@
 #include "Grid.h"
 namespace Grid {

-CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
+CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
 {
  _processors = processors;
  _ndimension = processors.size();
--- a/lib/cshift/Grid_cshift_common.h
+++ b/lib/cshift/Grid_cshift_common.h
--- a/lib/cshift/Grid_cshift_mpi.h
+++ b/lib/cshift/Grid_cshift_mpi.h
--- a/lib/cshift/Grid_cshift_none.h
+++ b/lib/cshift/Grid_cshift_none.h
--- a/lib/lattice/Grid_lattice_ET.h
+++ b/lib/lattice/Grid_lattice_ET.h
--- a/lib/lattice/Grid_lattice_arith.h
+++ b/lib/lattice/Grid_lattice_arith.h
--- a/lib/lattice/Grid_lattice_base.h
+++ b/lib/lattice/Grid_lattice_base.h
@@ -283,24 +283,23 @@ PARALLEL_FOR_LOOP



-#include <lattice/Grid_lattice_conformable.h>
+#include <lattice/Lattice_conformable.h>
 #define GRID_LATTICE_EXPRESSION_TEMPLATES
 #ifdef  GRID_LATTICE_EXPRESSION_TEMPLATES
-#include <lattice/Grid_lattice_ET.h>
+#include <lattice/Lattice_ET.h>
 #else 
-#include <lattice/Grid_lattice_overload.h>
+#include <lattice/Lattice_overload.h>
 #endif
-#include <lattice/Grid_lattice_arith.h>
-#include <lattice/Grid_lattice_trace.h>
-#include <lattice/Grid_lattice_transpose.h>
-#include <lattice/Grid_lattice_local.h>
-#include <lattice/Grid_lattice_reduction.h>
-#include <lattice/Grid_lattice_peekpoke.h>
-#include <lattice/Grid_lattice_reality.h>
-#include <Grid_extract.h>
-#include <lattice/Grid_lattice_coordinate.h>
-#include <lattice/Grid_lattice_rng.h>
-#include <lattice/Grid_lattice_transfer.h>
+#include <lattice/Lattice_arith.h>
+#include <lattice/Lattice_trace.h>
+#include <lattice/Lattice_transpose.h>
+#include <lattice/Lattice_local.h>
+#include <lattice/Lattice_reduction.h>
+#include <lattice/Lattice_peekpoke.h>
+#include <lattice/Lattice_reality.h>
+#include <lattice/Lattice_coordinate.h>
+#include <lattice/Lattice_rng.h>
+#include <lattice/Lattice_transfer.h>



--- a/lib/lattice/Grid_lattice_comparison.h
+++ b/lib/lattice/Grid_lattice_comparison.h
--- a/lib/lattice/Grid_lattice_conformable.h
+++ b/lib/lattice/Grid_lattice_conformable.h
--- a/lib/lattice/Grid_lattice_coordinate.h
+++ b/lib/lattice/Grid_lattice_coordinate.h
--- a/lib/lattice/Grid_lattice_local.h
+++ b/lib/lattice/Grid_lattice_local.h
--- a/lib/lattice/Grid_lattice_overload.h
+++ b/lib/lattice/Grid_lattice_overload.h
--- a/lib/lattice/Grid_lattice_peekpoke.h
+++ b/lib/lattice/Grid_lattice_peekpoke.h
--- a/lib/lattice/Grid_lattice_reality.h
+++ b/lib/lattice/Grid_lattice_reality.h
@@ -48,5 +48,16 @@ PARALLEL_FOR_LOOP
    }


+    template<class vobj> inline auto Ta(const Lattice<vobj> &z) -> Lattice<decltype(Ta(z._odata[0]))>
+    {
+      Lattice<decltype(Ta(z._odata[0]))> ret(z._grid);
+PARALLEL_FOR_LOOP
+        for(int ss=0;ss<z._grid->oSites();ss++){
+            ret._odata[ss] = Ta(z._odata[ss]);
+        }
+      return ret;
+    }
+
+
 }
 #endif
--- a/lib/lattice/Grid_lattice_reduction.h
+++ b/lib/lattice/Grid_lattice_reduction.h
--- a/lib/lattice/Grid_lattice_rng.h
+++ b/lib/lattice/Grid_lattice_rng.h
--- a/lib/lattice/Grid_lattice_trace.h
+++ b/lib/lattice/Grid_lattice_trace.h
--- a/lib/lattice/Grid_lattice_transfer.h
+++ b/lib/lattice/Grid_lattice_transfer.h
--- a/lib/lattice/Grid_lattice_transpose.h
+++ b/lib/lattice/Grid_lattice_transpose.h
--- a/lib/lattice/Grid_lattice_where.h
+++ b/lib/lattice/Grid_lattice_where.h
--- a/lib/math/Grid_math_arith.h
+++ b/lib/math/Grid_math_arith.h
@@ -1,11 +0,0 @@
-#ifndef GRID_MATH_ARITH_H
-#define GRID_MATH_ARITH_H
-
-#include <math/Grid_math_arith_add.h>
-#include <math/Grid_math_arith_sub.h>
-#include <math/Grid_math_arith_mac.h>
-#include <math/Grid_math_arith_mul.h>
-#include <math/Grid_math_arith_scalar.h>
-
-#endif
-
--- a/lib/parallelIO/GridNerscIO.h
+++ b/lib/parallelIO/GridNerscIO.h
--- a/lib/qcd/LinalgUtils.h
+++ b/lib/qcd/LinalgUtils.h
@@ -0,0 +1,113 @@
+#ifndef GRID_QCD_LINALG_UTILS_H
+#define GRID_QCD_LINALG_UTILS_H
+
+namespace Grid{
+namespace QCD{
+////////////////////////////////////////////////////////////////////////
+//This file brings additional linear combination assist that is helpful
+//to QCD such as chiral projectors and spin matrices applied to one of the inputs.
+//These routines support five-D chiral fermions and contain s-subslice indexing 
+//on the 5d (rb4d) checkerboarded lattices
+////////////////////////////////////////////////////////////////////////
+template<class vobj> 
+void axpby_ssp(Lattice<vobj> &z, RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp = a*x._odata[ss+s]+b*y._odata[ss+sp];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void ag5xpby_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    multGamma5(tmp(),a*x._odata[ss+s]());
+    tmp = tmp + b*y._odata[ss+sp];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void axpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    multGamma5(tmp(),b*y._odata[ss+sp]());
+    tmp = tmp + a*x._odata[ss+s];
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void ag5xpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp1;
+    vobj tmp2;
+    tmp1 = a*x._odata[ss+s]+b*y._odata[ss+sp];
+    multGamma5(tmp2(),tmp1());
+    vstream(z._odata[ss+s],tmp2);
+  }
+}
+
+template<class vobj> 
+void axpby_ssp_pminus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    spProj5m(tmp,y._odata[ss+sp]);
+    tmp = a*x._odata[ss+s]+b*tmp;
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+
+template<class vobj> 
+void axpby_ssp_pplus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,y);
+  conformable(x,z);
+  GridBase *grid=x._grid;
+  int Ls = grid->_rdimensions[0];
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
+    vobj tmp;
+    spProj5p(tmp,y._odata[ss+sp]);
+    tmp = a*x._odata[ss+s]+b*tmp;
+    vstream(z._odata[ss+s],tmp);
+  }
+}
+}}
+#endif 
--- a/lib/qcd/QCD.h
+++ b/lib/qcd/QCD.h
@@ -307,8 +307,10 @@ namespace QCD {
 }   //namespace QCD
 } // Grid

+#include <qcd/SpaceTimeGrid.h>
 #include <qcd/Dirac.h>
 #include <qcd/TwoSpinor.h>
+#include <qcd/LinalgUtils.h>
 #include <qcd/action/Actions.h>

 #endif
--- a/lib/qcd/SpaceTimeGrid.cc
+++ b/lib/qcd/SpaceTimeGrid.cc
@@ -0,0 +1,52 @@
+#include <Grid.h>
+
+namespace Grid { 
+  namespace QCD {
+
+/////////////////////////////////////////////////////////////////
+// Public interface
+/////////////////////////////////////////////////////////////////
+GridCartesian *SpaceTimeGrid::makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi)
+{
+  return new GridCartesian(latt,simd,mpi); 
+}
+GridRedBlackCartesian *SpaceTimeGrid::makeFourDimRedBlackGrid(const GridCartesian *FourDimGrid)
+{
+  return new GridRedBlackCartesian(FourDimGrid); 
+}
+
+GridCartesian         *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid)
+{
+  int N4=FourDimGrid->_ndimension;
+
+  std::vector<int> latt5(1,Ls);
+  std::vector<int> simd5(1,1);
+  std::vector<int>  mpi5(1,1);
+  
+  for(int d=0;d<N4;d++){
+    latt5.push_back(FourDimGrid->_fdimensions[d]);
+    simd5.push_back(FourDimGrid->_simd_layout[d]);
+     mpi5.push_back(FourDimGrid->_processors[d]);
+  }
+  return new GridCartesian(latt5,simd5,mpi5); 
+}
+
+GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid)
+{
+  int N4=FourDimGrid->_ndimension;
+  int cbd=1;
+  std::vector<int> latt5(1,Ls);
+  std::vector<int> simd5(1,1);
+  std::vector<int>  mpi5(1,1);
+  std::vector<int>   cb5(1,0);
+    
+  for(int d=0;d<N4;d++){
+    latt5.push_back(FourDimGrid->_fdimensions[d]);
+    simd5.push_back(FourDimGrid->_simd_layout[d]);
+     mpi5.push_back(FourDimGrid->_processors[d]);
+      cb5.push_back(  1);
+    }
+  return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd); 
+}
+
+}}
--- a/lib/qcd/SpaceTimeGrid.h
+++ b/lib/qcd/SpaceTimeGrid.h
@@ -0,0 +1,18 @@
+#ifndef GRID_QCD_SPACE_TIME_GRID_H
+#define GRID_QCD_SPACE_TIME_GRID_H
+namespace Grid {
+namespace QCD {
+
+class SpaceTimeGrid {
+ public:
+
+  static GridCartesian         *makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi);
+  static GridRedBlackCartesian *makeFourDimRedBlackGrid       (const GridCartesian *FourDimGrid);
+  static GridCartesian         *makeFiveDimGrid        (int Ls,const GridCartesian *FourDimGrid);
+  static GridRedBlackCartesian *makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid);
+
+};
+
+}}
+
+#endif
--- a/lib/qcd/action/Actions.h
+++ b/lib/qcd/action/Actions.h
@@ -1,10 +1,100 @@
 #ifndef GRID_QCD_ACTIONS_H
 #define GRID_QCD_ACTIONS_H

-#include <qcd/action/fermion/FermionAction.h>
-#include <qcd/action/fermion/WilsonCompressor.h>
-#include <qcd/action/fermion/WilsonKernels.h>
+
+// Some reorganisation likely required as both Chroma and IroIro
+// are separating the concept of the operator from that of action.
+//
+// The FermAction contains methods to create 
+//
+// * Linear operators             (Hermitian and non-hermitian)  .. my LinearOperator
+// * System solvers               (Hermitian and non-hermitian)  .. my OperatorFunction
+// * MultiShift System solvers    (Hermitian and non-hermitian)  .. my OperatorFunction
+
+
+////////////////////////////////////////////
+// Abstract base interface
+////////////////////////////////////////////
+#include <qcd/action/fermion/FermionOperator.h>
+
+////////////////////////////////////////////
+// Utility functions
+////////////////////////////////////////////
+#include <qcd/action/fermion/WilsonCompressor.h>     //used by all wilson type fermions
+#include <qcd/action/fermion/WilsonKernels.h>        //used by all wilson type fermions
+
+////////////////////////////////////////////
+// 4D formulations
+////////////////////////////////////////////
 #include <qcd/action/fermion/WilsonFermion.h>
-#include <qcd/action/fermion/FiveDimWilsonFermion.h>
+//#include <qcd/action/fermion/CloverFermion.h>
+
+////////////////////////////////////////////
+// 5D formulations...
+////////////////////////////////////////////
+
+#include <qcd/action/fermion/WilsonFermion5D.h> // used by all 5d overlap types
+
+//////////
+// Cayley
+//////////
+#include <qcd/action/fermion/CayleyFermion5D.h>
+
+#include <qcd/action/fermion/DomainWallFermion.h>
+#include <qcd/action/fermion/DomainWallFermion.h>
+
+#include <qcd/action/fermion/MobiusFermion.h>
+#include <qcd/action/fermion/ScaledShamirFermion.h>
+#include <qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h>
+
+#include <qcd/action/fermion/MobiusZolotarevFermion.h>
+#include <qcd/action/fermion/ShamirZolotarevFermion.h>
+#include <qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h>
+
+//////////////////////
+// Continued fraction
+//////////////////////
+#include <qcd/action/fermion/ContinuedFractionFermion5D.h>
+#include <qcd/action/fermion/OverlapWilsonContfracTanhFermion.h>
+#include <qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h>
+
+//////////////////////
+// Partial fraction
+//////////////////////
+#include <qcd/action/fermion/PartialFractionFermion5D.h>
+#include <qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h>
+#include <qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h>
+
+
+    // Chroma interface defining FermionAction
+    /*
+     template<typename T, typename P, typename Q>  class FermAct4D : public FermionAction<T,P,Q>
+     virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
+     virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
+     virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
+     virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
+     virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
+     virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
+     class DiffFermAct4D : public FermAct4D<T,P,Q>
+     virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
+     virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
+    */
+
+
+    // Chroma interface defining GaugeAction
+    /*
+      template<typename P, typename Q>   class GaugeAction
+  virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
+  virtual GaugeState<P,Q>* createState(const Q& q) const
+  virtual const GaugeBC<P,Q>& getGaugeBC() const
+  virtual const Set& getSet(void) const = 0;
+  virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const 
+  virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
+
+  class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
+  typedef multi1d<LatticeColorMatrix>  P;
+    */

 #endif
--- a/lib/qcd/action/fermion/.dirstamp
+++ b/lib/qcd/action/fermion/.dirstamp
--- a/lib/qcd/action/fermion/CayleyFermion5D.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5D.cc
@@ -0,0 +1,346 @@
+#include <Grid.h>
+namespace Grid {
+namespace QCD {
+
+ CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu,
+				  GridCartesian         &FiveDimGrid,
+				  GridRedBlackCartesian &FiveDimRedBlackGrid,
+				  GridCartesian         &FourDimGrid,
+				  GridRedBlackCartesian &FourDimRedBlackGrid,
+				  RealD _mass,RealD _M5) :
+   WilsonFermion5D(_Umu,
+		   FiveDimGrid,
+		   FiveDimRedBlackGrid,
+		   FourDimGrid,
+		   FourDimRedBlackGrid,_M5),
+   mass(_mass)
+ {
+ }
+
+  // override multiply
+  RealD CayleyFermion5D::M    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion Din(psi._grid);
+
+    // Assemble Din
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	Din = bs psi[s] + cs[s] psi[s+1}
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	//      Din+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pplus (Din,1.0,Din,-mass*cs[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(Din,bs[s],psi,-mass*cs[s],psi,s,0);
+	axpby_ssp_pplus (Din,1.0,Din,cs[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
+	axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
+      }
+    }
+
+    DW(Din,chi,DaggerNo);
+    // ((b D_W + D_w hop terms +1) on s-diag
+    axpby(chi,1.0,1.0,chi,psi); 
+
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ){
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,Ls-1);
+      } else if ( s==(Ls-1)) {
+	axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,0);
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
+      }
+    }
+    return norm2(chi);
+  }
+
+  RealD CayleyFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Under adjoint
+    //D1+        D1- P-    ->   D1+^dag   P+ D2-^dag
+    //D2- P+     D2+            P-D1-^dag D2+dag
+
+    LatticeFermion Din(psi._grid);
+    // Apply Dw
+    DW(psi,Din,DaggerYes); 
+
+    for(int s=0;s<Ls;s++){
+      // Collect the terms in DW
+      //	Chi = bs Din[s] + cs[s] Din[s+1}
+      //    Chi+= -mass*cs[s] psi[s+1}
+      if ( s==0 ) {
+	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,-mass*cs[Ls-1],Din,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus (chi,bs[s],Din,-mass*cs[0],Din,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
+      } else {
+	axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
+      }
+      // Collect the terms indept of DW
+      if ( s==0 ){
+	axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,Ls-1);
+      } else if ( s==(Ls-1)) {
+	axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,1.0,chi,-1.0,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
+      }
+    }
+    // ((b D_W + D_w hop terms +1) on s-diag
+    axpby (chi,1.0,1.0,chi,psi); 
+    return norm2(chi);
+  }
+
+  // half checkerboard operations
+  void CayleyFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion tmp(psi._grid);
+    // Assemble the 5d matrix
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	//	tmp = bs psi[s] + cs[s] psi[s+1}
+	//      tmp+= -mass*cs[s] psi[s+1}
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi ,s, s+1);
+	axpby_ssp_pplus(tmp,1.0,tmp,mass*ceo[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(tmp,beo[s],psi,mass*ceo[s],psi,s,0);
+	axpby_ssp_pplus(tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi,s,s+1);
+	axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
+      }
+    }
+    // Apply 4d dslash
+    if ( psi.checkerboard == Odd ) {
+      DhopEO(tmp,chi,DaggerNo);
+    } else {
+      DhopOE(tmp,chi,DaggerNo);
+    }
+  }
+
+  void CayleyFermion5D::MeooeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    LatticeFermion tmp(psi._grid);
+    // Apply 4d dslash
+    if ( psi.checkerboard == Odd ) {
+      DhopEO(psi,tmp,DaggerYes);
+    } else {
+      DhopOE(psi,tmp,DaggerYes);
+    }
+    // Assemble the 5d matrix
+    for(int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	axpby_ssp_pplus(chi,beo[s],tmp,   -ceo[s+1]  ,tmp,s,s+1);
+	axpby_ssp_pminus(chi,   1.0,chi,mass*ceo[Ls-1],tmp,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus(chi,beo[s],tmp,mass*ceo[0],tmp,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-ceo[s-1],tmp,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,beo[s],tmp,-ceo[s+1],tmp,s,s+1);
+	axpby_ssp_pminus(chi,1.0   ,chi,-ceo[s-1],tmp,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    for (int s=0;s<Ls;s++){
+      if ( s==0 ) {
+	axpby_ssp_pminus(chi,bee[s],psi ,-cee[s],psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,mass*cee[s],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pminus(chi,bee[s],psi,mass*cee[s],psi,s,0);
+	axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
+      } else {
+	axpby_ssp_pminus(chi,bee[s],psi,-cee[s],psi,s,s+1);
+	axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::MooeeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    for (int s=0;s<Ls;s++){
+      // Assemble the 5d matrix
+      if ( s==0 ) {
+	axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1]  ,psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0,chi,mass*cee[Ls-1],psi,s,Ls-1);
+      } else if ( s==(Ls-1)) { 
+	axpby_ssp_pplus(chi,bee[s],psi,mass*cee[0],psi,s,0);
+	axpby_ssp_pminus(chi,1.0,chi,-cee[s-1],psi,s,s-1);
+      } else {
+	axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1],psi,s,s+1);
+	axpby_ssp_pminus(chi,1.0   ,chi,-cee[s-1],psi,s,s-1);
+      }
+    }
+  }
+
+  void CayleyFermion5D::MooeeInv    (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Apply (L^{\prime})^{-1}
+    axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+    for (int s=1;s<Ls;s++){
+      axpby_ssp_pplus(chi,1.0,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
+    }
+    // L_m^{-1} 
+    for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
+      axpby_ssp_pminus(chi,1.0,chi,-leem[s],chi,Ls-1,s);
+    }
+    // U_m^{-1} D^{-1}
+    for (int s=0;s<Ls-1;s++){
+      // Chi[s] + 1/d chi[s] 
+      axpby_ssp_pplus(chi,1.0/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
+    }	
+    axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+    
+    // Apply U^{-1}
+    for (int s=Ls-2;s>=0;s--){
+      axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1);  // chi[Ls]
+    }
+  }
+
+  void CayleyFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
+  {
+    // Apply (U^{\prime})^{-dagger}
+    axpby_ssp (chi,1.0,psi,     0.0,psi,0,0);      // chi[0]=psi[0]
+    for (int s=1;s<Ls;s++){
+      axpby_ssp_pminus(chi,1.0,psi,-uee[s-1],chi,s,s-1);
+    }
+    // U_m^{-\dagger} 
+    for (int s=0;s<Ls-1;s++){
+      axpby_ssp_pplus(chi,1.0,chi,-ueem[s],chi,Ls-1,s);
+    }
+    // L_m^{-\dagger} D^{-dagger}
+    for (int s=0;s<Ls-1;s++){
+      axpby_ssp_pminus(chi,1.0/dee[s],chi,-leem[s]/dee[Ls-1],chi,s,Ls-1);
+    }	
+    axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable 
+    
+    // Apply L^{-dagger}
+    for (int s=Ls-2;s>=0;s--){
+      axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1);  // chi[Ls]
+    }
+  }
+  
+  // Tanh
+  void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c)
+  {
+    SetCoefficientsZolotarev(1.0,zdata,b,c);
+
+  }
+  //Zolo
+  void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c)
+  {
+
+    ///////////////////////////////////////////////////////////
+    // The Cayley coeffs (unprec)
+    ///////////////////////////////////////////////////////////
+    omega.resize(Ls);
+    bs.resize(Ls);
+    cs.resize(Ls);
+    as.resize(Ls);
+    
+    // 
+    // Ts = (    [bs+cs]Dw        )^-1 (    (bs+cs) Dw         )
+    //     -(g5  -------       -1 )    ( g5 ---------     + 1  )
+    //      (   {2+(bs-cs)Dw}     )    (    2+(bs-cs) Dw       )
+    //
+    //  bs = 1/2( (1/omega_s + 1)*b + (1/omega - 1)*c ) = 1/2(  1/omega(b+c) + (b-c) )
+    //  cs = 1/2( (1/omega_s - 1)*b + (1/omega + 1)*c ) = 1/2(  1/omega(b+c) - (b-c) )
+    //
+    // bs+cs = 0.5*( 1/omega(b+c) + (b-c) + 1/omega(b+c) - (b-c) ) = 1/omega(b+c)
+    // bs-cs = 0.5*( 1/omega(b+c) + (b-c) - 1/omega(b+c) + (b-c) ) = b-c
+    //
+    // So 
+    //
+    // Ts = (    [b+c]Dw/omega_s    )^-1 (    (b+c) Dw /omega_s        )
+    //     -(g5  -------         -1 )    ( g5 ---------           + 1  )
+    //      (   {2+(b-c)Dw}         )    (    2+(b-c) Dw               )
+    //
+    // Ts = (    [b+c]Dw            )^-1 (    (b+c) Dw                 )
+    //     -(g5  -------    -omega_s)    ( g5 ---------      + omega_s )
+    //      (   {2+(b-c)Dw}         )    (    2+(b-c) Dw               )
+    // 
+    
+    double bpc = b+c;
+    double bmc = b-c;
+    for(int i=0; i < Ls; i++){
+      as[i] = 1.0;
+      omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code
+      bs[i] = 0.5*(bpc/omega[i] + bmc);
+      cs[i] = 0.5*(bpc/omega[i] - bmc);
+    }
+
+    ////////////////////////////////////////////////////////
+    // Constants for the preconditioned matrix Cayley form
+    ////////////////////////////////////////////////////////
+    bee.resize(Ls);
+    cee.resize(Ls);
+    beo.resize(Ls);
+    ceo.resize(Ls);
+    
+    for(int i=0;i<Ls;i++){
+      bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
+      cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
+      beo[i]=as[i]*bs[i];
+      ceo[i]=-as[i]*cs[i];
+    }
+
+    aee.resize(Ls);
+    aeo.resize(Ls);
+    for(int i=0;i<Ls;i++){
+      aee[i]=cee[i];
+      aeo[i]=ceo[i];
+    }
+
+    //////////////////////////////////////////
+    // LDU decomposition of eeoo
+    //////////////////////////////////////////
+    dee.resize(Ls);
+    lee.resize(Ls);
+    leem.resize(Ls);
+    uee.resize(Ls);
+    ueem.resize(Ls);
+    
+    for(int i=0;i<Ls;i++){
+      
+      dee[i] = bee[i];
+      
+      if ( i < Ls-1 ) {
+	
+	lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
+	    
+	leem[i]=mass*cee[Ls-1]/bee[0];
+	for(int j=0;j<i;j++)  leem[i]*= aee[j]/bee[j+1];
+	
+	uee[i] =-aee[i]/bee[i];   // up-diag entry on the ith row
+	
+	ueem[i]=mass;
+	for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
+	ueem[i]*= aee[0]/bee[0];
+	    
+      } else { 
+	lee[i] =0.0;
+	leem[i]=0.0;
+	uee[i] =0.0;
+	ueem[i]=0.0;
+      }
+    }
+	
+    { 
+      double delta_d=mass*cee[Ls-1];
+      for(int j=0;j<Ls-1;j++) delta_d *= cee[j]/bee[j];
+      dee[Ls-1] += delta_d;
+    }
+  }
+
+}}
+
+
--- a/lib/qcd/action/fermion/CayleyFermion5D.h
+++ b/lib/qcd/action/fermion/CayleyFermion5D.h
@@ -0,0 +1,62 @@
+#ifndef  GRID_QCD_CAYLEY_FERMION_H
+#define  GRID_QCD_CAYLEY_FERMION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class CayleyFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operations
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Instantiatable(void)=0;
+      //    protected:
+      RealD mass;
+
+      // Cayley form Moebius (tanh and zolotarev)
+      std::vector<RealD> omega; 
+      std::vector<RealD> bs;    // S dependent coeffs
+      std::vector<RealD> cs;    
+      std::vector<RealD> as;    
+      // For preconditioning Cayley form
+      std::vector<RealD> bee;    
+      std::vector<RealD> cee;    
+      std::vector<RealD> aee;    
+      std::vector<RealD> beo;    
+      std::vector<RealD> ceo;    
+      std::vector<RealD> aeo;    
+      // LDU factorisation of the eeoo matrix
+      std::vector<RealD> lee;    
+      std::vector<RealD> leem;    
+      std::vector<RealD> uee;    
+      std::vector<RealD> ueem;    
+      std::vector<RealD> dee;    
+
+      // Constructors
+      CayleyFermion5D(LatticeGaugeField &_Umu,
+		      GridCartesian         &FiveDimGrid,
+		      GridRedBlackCartesian &FiveDimRedBlackGrid,
+		      GridCartesian         &FourDimGrid,
+		      GridRedBlackCartesian &FourDimRedBlackGrid,
+		      RealD _mass,RealD _M5);
+
+    protected:
+      void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c);
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c);
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
+++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
@@ -0,0 +1,179 @@
+#include <Grid.h>
+
+namespace Grid {
+  namespace QCD {
+
+    void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale)
+    {
+      SetCoefficientsZolotarev(1.0/scale,zdata);
+    }
+    void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
+    {
+      // How to check Ls matches??
+      //      std::cout << Ls << " Ls"<<std::endl;
+      //      std::cout << zdata->n  << " - n"<<std::endl;
+      //      std::cout << zdata->da << " -da "<<std::endl;
+      //      std::cout << zdata->db << " -db"<<std::endl;
+      //      std::cout << zdata->dn << " -dn"<<std::endl;
+      //      std::cout << zdata->dd << " -dd"<<std::endl;
+
+      assert(zdata->db==Ls);// Beta has Ls coeffs
+
+      R=(1+this->mass)/(1-this->mass);
+
+      Beta.resize(Ls);
+      cc.resize(Ls);
+      cc_d.resize(Ls);
+      sqrt_cc.resize(Ls);
+      for(int i=0; i < Ls ; i++){
+	Beta[i] = zdata -> beta[i];
+	cc[i] = 1.0/Beta[i];
+	cc_d[i]=sqrt(cc[i]);
+      }
+    
+      cc_d[Ls-1]=1.0;
+      for(int i=0; i < Ls-1 ; i++){
+	sqrt_cc[i]= sqrt(cc[i]*cc[i+1]);
+      }    
+      sqrt_cc[Ls-2]=sqrt(cc[Ls-2]);
+
+
+      ZoloHiInv =1.0/zolo_hi;
+      dw_diag = (4.0-M5)*ZoloHiInv;
+    
+      See.resize(Ls);
+      Aee.resize(Ls);
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	Aee[s] = sign * Beta[s] * dw_diag;
+	sign   = - sign;
+      }
+      Aee[Ls-1] += R;
+    
+      See[0] = Aee[0];
+      for(int s=1;s<Ls;s++){
+	See[s] = Aee[s] - 1.0/See[s-1];
+      }
+      for(int s=0;s<Ls;s++){
+	std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
+      }
+    }
+
+
+
+    RealD  ContinuedFractionFermion5D::M           (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      LatticeFermion D(psi._grid);
+
+      DW(psi,D,DaggerNo); 
+
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==0 ) {
+	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*ZoloHiInv,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
+	} else if ( s==(Ls-1) ){
+	  RealD R=(1.0+mass)/(1.0-mass);
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,D,sqrt_cc[s-1],psi,s,s-1);
+	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,D,sqrt_cc[s],psi,s,s+1);
+  	  axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
+	}
+	sign=-sign; 
+      }
+      return norm2(chi);
+    }
+    RealD  ContinuedFractionFermion5D::Mdag        (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag
+      // The rest of matrix is symmetric.
+      // Can ignore "dag"
+      return M(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::Meooe       (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // Apply 4d dslash
+      if ( psi.checkerboard == Odd ) {
+	DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      } else {
+	DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
+      }
+      
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==(Ls-1) ){
+	  ag5xpby_ssp(chi,Beta[s]*ZoloHiInv,chi,0.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*ZoloHiInv,chi,0.0,chi,s,s);
+	}
+	sign=-sign; 
+      }
+    }
+    void   ContinuedFractionFermion5D::MeooeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      Meooe(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      int sign=1;
+      for(int s=0;s<Ls;s++){
+	if ( s==0 ) {
+	  ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*dw_diag,psi,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
+	} else if ( s==(Ls-1) ){
+	  // Drop the CC here.
+	  double R=(1+mass)/(1-mass);
+	  ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1);
+	  ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
+	} else {
+	  ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*dw_diag,psi,sqrt_cc[s],psi,s,s+1);
+	  axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
+	}
+	sign=-sign; 
+      }
+    }
+
+    void   ContinuedFractionFermion5D::MooeeDag    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      Mooee(psi,chi);
+    }
+    void   ContinuedFractionFermion5D::MooeeInv    (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      // Apply Linv
+      axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); 
+      for(int s=1;s<Ls;s++){
+	axpbg5y_ssp(chi,1.0/cc_d[s],psi,-1.0/See[s-1],chi,s,s-1);
+      }
+      // Apply Dinv
+      for(int s=0;s<Ls;s++){
+	ag5xpby_ssp(chi,1.0/See[s],chi,0.0,chi,s,s); //only appearance of See[0]
+      }
+      // Apply Uinv = (Linv)^T
+      axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,Ls-1,Ls-1);
+      for(int s=Ls-2;s>=0;s--){
+	axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1);
+      }
+    }
+    void   ContinuedFractionFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
+    {
+      MooeeInv(psi,chi);
+    }
+    
+    // Constructors
+    ContinuedFractionFermion5D::ContinuedFractionFermion5D(
+							   LatticeGaugeField &_Umu,
+							   GridCartesian         &FiveDimGrid,
+							   GridRedBlackCartesian &FiveDimRedBlackGrid,
+							   GridCartesian         &FourDimGrid,
+							   GridRedBlackCartesian &FourDimRedBlackGrid,
+							   RealD _mass,RealD M5) :
+      WilsonFermion5D(_Umu,
+		      FiveDimGrid, FiveDimRedBlackGrid,
+		      FourDimGrid, FourDimRedBlackGrid,M5),
+      mass(_mass)
+    {
+      assert((Ls&0x1)==1); // Odd Ls required
+    }
+
+  }
+}
+
--- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h
+++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h
@@ -0,0 +1,58 @@
+#ifndef  GRID_QCD_CONTINUED_FRACTION_H
+#define  GRID_QCD_CONTINUED_FRACTION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ContinuedFractionFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operaions
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+
+      //      virtual void   Instantiatable(void)=0;
+      virtual void   Instantiatable(void) =0;
+
+      // Constructors
+      ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
+				 GridCartesian         &FiveDimGrid,
+				 GridRedBlackCartesian &FiveDimRedBlackGrid,
+				 GridCartesian         &FourDimGrid,
+				 GridRedBlackCartesian &FourDimRedBlackGrid,
+				 RealD _mass,RealD M5);
+
+    protected:
+
+      void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
+      void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);;
+
+      // Cont frac
+      RealD dw_diag;
+      RealD mass;
+      RealD R;
+      RealD ZoloHiInv;
+      std::vector<double> Beta;
+      std::vector<double> cc;;
+      std::vector<double> cc_d;;
+      std::vector<double> sqrt_cc;
+      std::vector<double> See;
+      std::vector<double> Aee;
+
+    };
+
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/DomainWallFermion.h
+++ b/lib/qcd/action/fermion/DomainWallFermion.h
@@ -0,0 +1,47 @@
+#ifndef  GRID_QCD_DOMAIN_WALL_FERMION_H
+#define  GRID_QCD_DOMAIN_WALL_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class DomainWallFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+      DomainWallFermion(LatticeGaugeField &_Umu,
+			GridCartesian         &FiveDimGrid,
+			GridRedBlackCartesian &FiveDimRedBlackGrid,
+			GridCartesian         &FourDimGrid,
+			GridRedBlackCartesian &FourDimRedBlackGrid,
+			RealD _mass,RealD _M5) : 
+
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = 1.0;
+
+	Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
+	assert(zdata->n==this->Ls);
+	
+	std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);
+
+	Approx::zolotarev_free(zdata);
+      }
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/FermionOperator.h
+++ b/lib/qcd/action/fermion/FermionOperator.h
@@ -1,5 +1,5 @@
-#ifndef  GRID_QCD_WILSON_DOP_H
-#define  GRID_QCD_WILSON_DOP_H
+#ifndef  GRID_QCD_FERMION_OPERATOR_H
+#define  GRID_QCD_FERMION_OPERATOR_H

 namespace Grid {

@@ -11,7 +11,7 @@ namespace Grid {
    // Think about multiple representations
    //////////////////////////////////////////////////////////////////////////////
    template<class FermionField,class GaugeField>
-    class FermionAction : public CheckerBoardedSparseMatrixBase<FermionField>
+    class FermionOperator : public CheckerBoardedSparseMatrixBase<FermionField>
    {
    public:

@@ -40,6 +40,7 @@ namespace Grid {
      virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
      virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;

+
    };

  }
--- a/lib/qcd/action/fermion/MobiusFermion.h
+++ b/lib/qcd/action/fermion/MobiusFermion.h
@@ -0,0 +1,49 @@
+#ifndef  GRID_QCD_MOBIUS_FERMION_H
+#define  GRID_QCD_MOBIUS_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class MobiusFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+      MobiusFermion(LatticeGaugeField &_Umu,
+		    GridCartesian         &FiveDimGrid,
+		    GridRedBlackCartesian &FiveDimRedBlackGrid,
+		    GridCartesian         &FourDimGrid,
+		    GridRedBlackCartesian &FourDimRedBlackGrid,
+		    RealD _mass,RealD _M5,
+		    RealD b, RealD c) : 
+      
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = 1.0;
+
+	std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
+	Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
+	assert(zdata->n==this->Ls);
+	
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c);
+
+	Approx::zolotarev_free(zdata);
+ 
+      }
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/MobiusZolotarevFermion.h
+++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h
@@ -0,0 +1,50 @@
+#ifndef  GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
+#define  GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class MobiusZolotarevFermion : public CayleyFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+       MobiusZolotarevFermion(LatticeGaugeField &_Umu,
+			      GridCartesian         &FiveDimGrid,
+			      GridRedBlackCartesian &FiveDimRedBlackGrid,
+			      GridCartesian         &FourDimGrid,
+			      GridRedBlackCartesian &FourDimRedBlackGrid,
+			      RealD _mass,RealD _M5,
+			      RealD b, RealD c,
+			      RealD lo, RealD hi) : 
+      
+      CayleyFermion5D(_Umu,
+		      FiveDimGrid,
+		      FiveDimRedBlackGrid,
+		      FourDimGrid,
+		      FourDimRedBlackGrid,_mass,_M5)
+
+      {
+	RealD eps = lo/hi;
+
+	Approx::zolotarev_data *zdata = Approx::zolotarev(eps,this->Ls,0);
+	assert(zdata->n==this->Ls);
+
+	std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
+	
+	// Call base setter
+	this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);
+ 
+	Approx::zolotarev_free(zdata);
+      }
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h
@@ -0,0 +1,34 @@
+#ifndef OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
+#define OVERLAP_WILSON_CAYLEY_TANH_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonCayleyTanhFermion : public MobiusFermion
+    {
+    public:
+
+      // Constructors
+    OverlapWilsonCayleyTanhFermion(LatticeGaugeField &_Umu,
+				   GridCartesian         &FiveDimGrid,
+				   GridRedBlackCartesian &FiveDimRedBlackGrid,
+				   GridCartesian         &FourDimGrid,
+				   GridRedBlackCartesian &FourDimRedBlackGrid,
+				   RealD _mass,RealD _M5,
+				   RealD scale) :
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      MobiusFermion(_Umu,
+		    FiveDimGrid,
+		    FiveDimRedBlackGrid,
+		    FourDimGrid,
+		    FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale)
+	{
+	}
+    };
+  }
+}
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h
@@ -0,0 +1,37 @@
+#ifndef  OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
+#define  OVERLAP_WILSON_CAYLEY_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion
+    {
+    public:
+
+      // Constructors
+
+    OverlapWilsonCayleyZolotarevFermion(LatticeGaugeField &_Umu,
+					GridCartesian         &FiveDimGrid,
+					GridRedBlackCartesian &FiveDimRedBlackGrid,
+					GridCartesian         &FourDimGrid,
+					GridRedBlackCartesian &FourDimRedBlackGrid,
+					RealD _mass,RealD _M5,
+					RealD lo, RealD hi) : 
+      // b+c=1.0, b-c = 0 <=> b =c = 1/2
+      MobiusZolotarevFermion(_Umu,
+			     FiveDimGrid,
+			     FiveDimRedBlackGrid,
+			     FourDimGrid,
+			     FourDimRedBlackGrid,_mass,_M5,0.5,0.5,lo,hi)
+
+      {}
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h
@@ -0,0 +1,40 @@
+#ifndef OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
+#define OVERLAP_WILSON_CONTFRAC_TANH_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonContFracTanhFermion(LatticeGaugeField &_Umu,
+				     GridCartesian         &FiveDimGrid,
+				     GridRedBlackCartesian &FiveDimRedBlackGrid,
+				     GridCartesian         &FourDimGrid,
+				     GridRedBlackCartesian &FourDimRedBlackGrid,
+				     RealD _mass,RealD _M5,
+				     RealD scale) :
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      ContinuedFractionFermion5D(_Umu,
+				 FiveDimGrid,
+				 FiveDimRedBlackGrid,
+				 FourDimGrid,
+				 FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+	  int nrational=Ls-1;// Even rational order
+	  Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham
+	  SetCoefficientsTanh(zdata,scale);
+	  Approx::zolotarev_free(zdata);
+	}
+    };
+  }
+}
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h
@@ -0,0 +1,44 @@
+#ifndef OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
+#define OVERLAP_WILSON_CONTFRAC_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonContFracZolotarevFermion(LatticeGaugeField &_Umu,
+					  GridCartesian         &FiveDimGrid,
+					  GridRedBlackCartesian &FiveDimRedBlackGrid,
+					  GridCartesian         &FourDimGrid,
+					  GridRedBlackCartesian &FourDimRedBlackGrid,
+					  RealD _mass,RealD _M5,
+					  RealD lo,RealD hi):
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      ContinuedFractionFermion5D(_Umu,
+				 FiveDimGrid,
+				 FiveDimRedBlackGrid,
+				 FourDimGrid,
+				 FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+
+	  int nrational=Ls;// Odd rational order
+	  RealD eps = lo/hi;
+
+	  Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0);
+	  SetCoefficientsZolotarev(hi,zdata);
+	  Approx::zolotarev_free(zdata);
+
+	}
+    };
+  }
+}
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h
@@ -0,0 +1,40 @@
+#ifndef OVERLAP_WILSON_PARTFRAC_TANH_FERMION_H
+#define OVERLAP_WILSON_PARTFRAC_TANH_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonPartialFractionTanhFermion(LatticeGaugeField &_Umu,
+					    GridCartesian         &FiveDimGrid,
+					    GridRedBlackCartesian &FiveDimRedBlackGrid,
+					    GridCartesian         &FourDimGrid,
+					    GridRedBlackCartesian &FourDimRedBlackGrid,
+					    RealD _mass,RealD _M5,
+					    RealD scale) :
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      PartialFractionFermion5D(_Umu,
+			       FiveDimGrid,
+			       FiveDimRedBlackGrid,
+			       FourDimGrid,
+			       FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+	  int nrational=Ls-1;// Even rational order
+	  Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham
+	  SetCoefficientsTanh(zdata,scale);
+	  Approx::zolotarev_free(zdata);
+	}
+    };
+  }
+}
+#endif
--- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h
+++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h
@@ -0,0 +1,44 @@
+#ifndef OVERLAP_WILSON_PARTFRAC_ZOLOTAREV_FERMION_H
+#define OVERLAP_WILSON_PARTFRAC_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D
+    {
+    public:
+
+      virtual void   Instantiatable(void){};
+      // Constructors
+    OverlapWilsonPartialFractionZolotarevFermion(LatticeGaugeField &_Umu,
+					  GridCartesian         &FiveDimGrid,
+					  GridRedBlackCartesian &FiveDimRedBlackGrid,
+					  GridCartesian         &FourDimGrid,
+					  GridRedBlackCartesian &FourDimRedBlackGrid,
+					  RealD _mass,RealD _M5,
+					  RealD lo,RealD hi):
+      
+      // b+c=scale, b-c = 0 <=> b =c = scale/2
+      PartialFractionFermion5D(_Umu,
+			       FiveDimGrid,
+			       FiveDimRedBlackGrid,
+			       FourDimGrid,
+			       FourDimRedBlackGrid,_mass,_M5)
+	{
+	  assert((Ls&0x1)==1); // Odd Ls required
+
+	  int nrational=Ls;// Odd rational order
+	  RealD eps = lo/hi;
+
+	  Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0);
+	  SetCoefficientsZolotarev(hi,zdata);
+	  Approx::zolotarev_free(zdata);
+
+	}
+    };
+  }
+}
+#endif
--- a/lib/qcd/action/fermion/PartialFractionFermion5D.cc
+++ b/lib/qcd/action/fermion/PartialFractionFermion5D.cc
@@ -0,0 +1,310 @@
+#include <Grid.h>
+namespace Grid {
+  namespace QCD {
+
+    void   PartialFractionFermion5D::Meooe_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
+    {
+      // this does both dag and undag but is trivial; make a common helper routing
+      int sign = dag ? (-1) : 1;
+
+      if ( psi.checkerboard == Odd ) {
+	DhopEO(psi,chi,DaggerNo);
+      } else {
+	DhopOE(psi,chi,DaggerNo);
+      }
+
+      int nblock=(Ls-1)/2;
+      for(int b=0;b<nblock;b++){
+	int s = 2*b;
+	ag5xpby_ssp(chi,-scale,chi,0.0,chi,s,s); 
+	ag5xpby_ssp(chi, scale,chi,0.0,chi,s+1,s+1); 
+      }
+      ag5xpby_ssp(chi,p[nblock]*scale/amax,chi,0.0,chi,Ls-1,Ls-1);
+    }
+
+    void   PartialFractionFermion5D::Mooee_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
+    {
+      // again dag and undag are trivially related
+      int sign = dag ? (-1) : 1;
+      
+      int nblock=(Ls-1)/2;
+      for(int b=0;b<nblock;b++){
+	
+	int s = 2*b;
+	RealD pp = p[nblock-1-b];
+	RealD qq = q[nblock-1-b];
+	
+	// Do each 2x2 block aligned at s and multiplies Dw site diagonal by G5 so Hw
+	ag5xpby_ssp(chi,-dw_diag*scale,psi,amax*sqrt(qq)*scale,psi, s  ,s+1); 
+	ag5xpby_ssp(chi, dw_diag*scale,psi,amax*sqrt(qq)*scale,psi, s+1,s);
+	axpby_ssp  (chi, 1.0, chi,sqrt(amax*pp)*scale*sign,psi,s+1,Ls-1);
+      }
+      
+      {
+	RealD R=(1+mass)/(1-mass);
+	//R g5 psi[Ls-1] + p[0] H
+	ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale*dw_diag/amax,psi,Ls-1,Ls-1);
+	
+	for(int b=0;b<nblock;b++){
+	  int s = 2*b+1;
+	  RealD pp = p[nblock-1-b];
+	  axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
+	}
+      }
+    }
+
+    void   PartialFractionFermion5D::MooeeInv_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
+    {
+      int sign = dag ? (-1) : 1;
+
+      LatticeFermion tmp(psi._grid);
+      
+      ///////////////////////////////////////////////////////////////////////////////////////
+      //Linv
+      ///////////////////////////////////////////////////////////////////////////////////////
+      int nblock=(Ls-1)/2;
+
+      axpy(chi,0.0,psi,psi); // Identity piece
+      
+      for(int b=0;b<nblock;b++){
+	int s = 2*b;
+	RealD pp = p[nblock-1-b];
+	RealD qq = q[nblock-1-b];
+	RealD coeff1=sign*sqrt(amax*amax*amax*pp*qq) / ( dw_diag*dw_diag + amax*amax* qq);
+	RealD coeff2=sign*sqrt(amax*pp)*dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
+	axpby_ssp  (chi,1.0,chi,coeff1,psi,Ls-1,s);
+	axpbg5y_ssp(chi,1.0,chi,coeff2,psi,Ls-1,s+1);
+      }
+      
+      ///////////////////////////////////////////////////////////////////////////////////////
+      //Dinv (note D isn't really diagonal -- just diagonal enough that we can still invert)
+      // Compute Seeinv (coeff of gamma5)
+      ///////////////////////////////////////////////////////////////////////////////////////
+      RealD R=(1+mass)/(1-mass);
+      RealD Seeinv = R + p[nblock]*dw_diag/amax;
+      for(int b=0;b<nblock;b++){
+	Seeinv += p[nblock-1-b]*dw_diag/amax / ( dw_diag*dw_diag/amax/amax + q[nblock-1-b]);
+      }    
+      Seeinv = 1.0/Seeinv;
+      
+      for(int b=0;b<nblock;b++){
+	int s = 2*b;
+	RealD pp = p[nblock-1-b];
+	RealD qq = q[nblock-1-b];
+	RealD coeff1=dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
+	RealD coeff2=amax*sqrt(qq) / ( dw_diag*dw_diag + amax*amax* qq);
+	ag5xpby_ssp  (tmp,-coeff1,chi,coeff2,chi,s,s+1);
+	ag5xpby_ssp  (tmp, coeff1,chi,coeff2,chi,s+1,s);
+      }
+      ag5xpby_ssp  (tmp, Seeinv,chi,0.0,chi,Ls-1,Ls-1);
+      
+      ///////////////////////////////////////////////////////////////////////////////////////
+      // Uinv
+      ///////////////////////////////////////////////////////////////////////////////////////
+      for(int b=0;b<nblock;b++){
+	int s = 2*b;
+	RealD pp = p[nblock-1-b];
+	RealD qq = q[nblock-1-b];
+	RealD coeff1=-sign*sqrt(amax*amax*amax*pp*qq) / ( dw_diag*dw_diag + amax*amax* qq);
+	RealD coeff2=-sign*sqrt(amax*pp)*dw_diag / ( dw_diag*dw_diag + amax*amax* qq); // Implicit g5 here
+	axpby_ssp  (chi,1.0/scale,tmp,coeff1/scale,tmp,s,Ls-1);
+	axpbg5y_ssp(chi,1.0/scale,tmp,coeff2/scale,tmp,s+1,Ls-1);
+      }
+      axpby_ssp  (chi, 1.0/scale,tmp,0.0,tmp,Ls-1,Ls-1);
+    }
+
+    void   PartialFractionFermion5D::M_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag)
+    {
+      LatticeFermion D(psi._grid);
+  
+      int sign = dag ? (-1) : 1;
+
+      // For partial frac Hw case (b5=c5=1) chroma quirkily computes
+      //
+      // Conventions for partfrac appear to be a mess.
+      // Tony's Nara lectures have
+      //
+      // BlockDiag(  H/p_i  1             | 1       )    
+      //          (  1      p_i H / q_i^2 | 0       )  
+      //           ---------------------------------
+      //           ( -1      0                | R  +p0 H  )
+      //
+      //Chroma     ( -2H    2sqrt(q_i)    |   0         )
+      //           (2 sqrt(q_i)   2H      |  2 sqrt(p_i) )
+      //           ---------------------------------
+      //           ( 0     -2 sqrt(p_i)   |  2 R gamma_5 + p0 2H
+      //
+      // Edwards/Joo/Kennedy/Wenger
+      //
+      // Here, the "beta's" selected by chroma to scale the unphysical bulk constraint fields
+      // incorporate the approx scale factor. This is obtained by propagating the
+      // scale on "H" out to the off diagonal elements as follows:
+      //
+      // BlockDiag(  H/p_i  1             | 1       ) 
+      //          (  1      p_i H / q_i^2 | 0       )  
+      //           ---------------------------------
+      //          ( -1      0                | R  + p_0 H  )
+      //
+      // becomes:
+      // BlockDiag(  H/ sp_i  1               | 1             ) 
+      //          (  1      sp_i H / s^2q_i^2 | 0             )  
+      //           ---------------------------------
+      //           ( -1      0                | R + p_0/s H   )
+      //
+      //
+      // This is implemented in Chroma by
+      //           p0' = p0/approxMax
+      //           p_i' = p_i*approxMax
+      //           q_i' = q_i*approxMax*approxMax
+      //
+      // After the equivalence transform is applied the matrix becomes
+      // 
+      //Chroma     ( -2H    sqrt(q'_i)    |   0         )
+      //           (sqrt(q'_i)   2H       |   sqrt(p'_i) )
+      //           ---------------------------------
+      //           ( 0     -sqrt(p'_i)    |  2 R gamma_5 + p'0 2H
+      //
+      //     =     ( -2H    sqrt(q_i)amax    |   0              )
+      //           (sqrt(q_i)amax   2H       |   sqrt(p_i*amax) )
+      //           ---------------------------------
+      //           ( 0     -sqrt(p_i)*amax   |  2 R gamma_5 + p0/amax 2H
+      //
+
+      DW(psi,D,DaggerNo); 
+
+      int nblock=(Ls-1)/2;
+      for(int b=0;b<nblock;b++){
+	
+	int s = 2*b;
+	double pp = p[nblock-1-b];
+	double qq = q[nblock-1-b];
+	
+	// Do each 2x2 block aligned at s and
+	ag5xpby_ssp(chi,-1.0*scale,D,amax*sqrt(qq)*scale,psi, s  ,s+1); // Multiplies Dw by G5 so Hw
+	ag5xpby_ssp(chi, 1.0*scale,D,amax*sqrt(qq)*scale,psi, s+1,s);
+	
+	// Pick up last column
+	axpby_ssp  (chi, 1.0, chi,sqrt(amax*pp)*scale*sign,psi,s+1,Ls-1);
+      }
+	
+      {
+	double R=(1+this->mass)/(1-this->mass);
+	//R g5 psi[Ls] + p[0] H
+	ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1);
+	
+	for(int b=0;b<nblock;b++){
+	  int s = 2*b+1;
+	  double pp = p[nblock-1-b];
+	  axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
+	}
+      }
+
+    }
+
+    RealD  PartialFractionFermion5D::M    (const LatticeFermion &in, LatticeFermion &out)
+    {
+      M_internal(in,out,DaggerNo);
+      return norm2(out);
+    }
+    RealD  PartialFractionFermion5D::Mdag (const LatticeFermion &in, LatticeFermion &out)
+    {
+      M_internal(in,out,DaggerYes);
+      return norm2(out);
+    }
+
+    void PartialFractionFermion5D::Meooe       (const LatticeFermion &in, LatticeFermion &out)
+    {
+      Meooe_internal(in,out,DaggerNo);
+    }
+    void PartialFractionFermion5D::MeooeDag    (const LatticeFermion &in, LatticeFermion &out)
+    {
+      Meooe_internal(in,out,DaggerYes);
+    }
+    void PartialFractionFermion5D::Mooee       (const LatticeFermion &in, LatticeFermion &out)
+    {
+      Mooee_internal(in,out,DaggerNo);
+    }
+    void PartialFractionFermion5D::MooeeDag    (const LatticeFermion &in, LatticeFermion &out)
+    {
+      Mooee_internal(in,out,DaggerYes);
+    }
+
+    void PartialFractionFermion5D::MooeeInv    (const LatticeFermion &in, LatticeFermion &out)
+    {
+      MooeeInv_internal(in,out,DaggerNo);
+    }
+    void PartialFractionFermion5D::MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)
+    {
+      MooeeInv_internal(in,out,DaggerYes);
+    }
+
+    void  PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){
+      SetCoefficientsZolotarev(1.0/scale,zdata);
+    }
+    void  PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){
+
+      // check on degree matching
+      //      std::cout << Ls << " Ls"<<std::endl;
+      //      std::cout << zdata->n  << " - n"<<std::endl;
+      //      std::cout << zdata->da << " -da "<<std::endl;
+      //      std::cout << zdata->db << " -db"<<std::endl;
+      //      std::cout << zdata->dn << " -dn"<<std::endl;
+      //      std::cout << zdata->dd << " -dd"<<std::endl;
+      assert(Ls == (2*zdata->da -1) );
+
+      // Part frac
+      //      RealD R;
+      R=(1+mass)/(1-mass);
+      dw_diag = (4.0-M5);
+
+      //      std::vector<RealD> p; 
+      //      std::vector<RealD> q;
+      p.resize(zdata->da);
+      q.resize(zdata->dd);
+	
+      for(int n=0;n<zdata->da;n++){
+	p[n] = zdata -> alpha[n];
+      }
+      for(int n=0;n<zdata->dd;n++){
+	q[n] = -zdata -> ap[n];
+      }
+      
+      scale= part_frac_chroma_convention ? 2.0 : 1.0; // Chroma conventions annoy me
+
+      amax=zolo_hi;
+    }
+
+      // Constructors
+    PartialFractionFermion5D::PartialFractionFermion5D(LatticeGaugeField &_Umu,
+						       GridCartesian         &FiveDimGrid,
+						       GridRedBlackCartesian &FiveDimRedBlackGrid,
+						       GridCartesian         &FourDimGrid,
+						       GridRedBlackCartesian &FourDimRedBlackGrid,
+						       RealD _mass,RealD M5) :
+      WilsonFermion5D(_Umu,
+		      FiveDimGrid, FiveDimRedBlackGrid,
+		      FourDimGrid, FourDimRedBlackGrid,M5),
+      mass(_mass)
+
+    {
+      assert((Ls&0x1)==1); // Odd Ls required
+      int nrational=Ls-1;
+
+
+      Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);
+
+      // NB: chroma uses a cast to "float" for the zolotarev range(!?).
+      // this creates a real difference in the operator which I do not like but we can replicate here
+      // to demonstrate compatibility
+      //      RealD eps = (zolo_lo / zolo_hi);
+      //      zdata = bfm_zolotarev(eps,nrational,0);
+      
+      SetCoefficientsTanh(zdata,1.0);
+
+      Approx::zolotarev_free(zdata);
+
+    }
+ 
+ }
+}
+
--- a/lib/qcd/action/fermion/PartialFractionFermion5D.h
+++ b/lib/qcd/action/fermion/PartialFractionFermion5D.h
@@ -0,0 +1,61 @@
+#ifndef  GRID_QCD_PARTIAL_FRACTION_H
+#define  GRID_QCD_PARTIAL_FRACTION_H
+
+namespace Grid {
+
+  namespace QCD {
+
+    class PartialFractionFermion5D : public WilsonFermion5D
+    {
+    public:
+
+      const int part_frac_chroma_convention=1;
+
+      void   Meooe_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
+      void   Mooee_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
+      void   MooeeInv_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
+      void   M_internal(const LatticeFermion &in, LatticeFermion &out,int dag);
+
+      // override multiply
+      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
+      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+
+      // half checkerboard operaions
+      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
+      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+
+      virtual void   Instantiatable(void) =0; // ensure no make-eee
+
+      // Constructors
+      PartialFractionFermion5D(LatticeGaugeField &_Umu,
+				    GridCartesian         &FiveDimGrid,
+				    GridRedBlackCartesian &FiveDimRedBlackGrid,
+				    GridCartesian         &FourDimGrid,
+				    GridRedBlackCartesian &FourDimRedBlackGrid,
+				    RealD _mass,RealD M5);
+
+    protected:
+
+      virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
+      virtual void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);
+
+      // Part frac
+      RealD mass;
+      RealD dw_diag;
+      RealD R;
+      RealD amax;
+      RealD scale;
+      std::vector<double> p; 
+      std::vector<double> q;
+
+    };
+
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/ScaledShamirFermion.h
+++ b/lib/qcd/action/fermion/ScaledShamirFermion.h
@@ -0,0 +1,37 @@
+#ifndef  GRID_QCD_SCALED_SHAMIR_FERMION_H
+#define  GRID_QCD_SCALED_SHAMIR_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ScaledShamirFermion : public MobiusFermion
+    {
+    public:
+
+      // Constructors
+    ScaledShamirFermion(LatticeGaugeField &_Umu,
+			GridCartesian         &FiveDimGrid,
+			GridRedBlackCartesian &FiveDimRedBlackGrid,
+			GridCartesian         &FourDimGrid,
+			GridRedBlackCartesian &FourDimRedBlackGrid,
+			RealD _mass,RealD _M5,
+			RealD scale) :
+      
+      // b+c=scale, b-c = 1 <=> 2b = scale+1; 2c = scale-1
+      MobiusFermion(_Umu,
+		    FiveDimGrid,
+		    FiveDimRedBlackGrid,
+		    FourDimGrid,
+		    FourDimRedBlackGrid,_mass,_M5,0.5*(scale+1.0),0.5*(scale-1.0))
+      {
+      }
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/ShamirZolotarevFermion.h
+++ b/lib/qcd/action/fermion/ShamirZolotarevFermion.h
@@ -0,0 +1,39 @@
+#ifndef  GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
+#define  GRID_QCD_SHAMIR_ZOLOTAREV_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    class ShamirZolotarevFermion : public MobiusZolotarevFermion
+    {
+    public:
+
+      // Constructors
+
+
+    ShamirZolotarevFermion(LatticeGaugeField &_Umu,
+			   GridCartesian         &FiveDimGrid,
+			   GridRedBlackCartesian &FiveDimRedBlackGrid,
+			   GridCartesian         &FourDimGrid,
+			   GridRedBlackCartesian &FourDimRedBlackGrid,
+			   RealD _mass,RealD _M5,
+			   RealD lo, RealD hi) : 
+      
+      // b+c = 1; b-c = 1 => b=1, c=0
+      MobiusZolotarevFermion(_Umu,
+			     FiveDimGrid,
+			     FiveDimRedBlackGrid,
+			     FourDimGrid,
+			     FourDimRedBlackGrid,_mass,_M5,1.0,0.0,lo,hi)
+      
+      {}
+
+    };
+
+  }
+}
+
+#endif
--- a/lib/qcd/action/fermion/WilsonFermion.cc
+++ b/lib/qcd/action/fermion/WilsonFermion.cc
@@ -9,9 +9,9 @@ const std::vector<int> WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
 int WilsonFermion::HandOptDslash;

 WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu,
-			   GridCartesian         &Fgrid,
-			   GridRedBlackCartesian &Hgrid, 
-			   double _mass) :
+			     GridCartesian         &Fgrid,
+			     GridRedBlackCartesian &Hgrid, 
+			     RealD _mass) :
  _grid(&Fgrid),
  _cbgrid(&Hgrid),
  Stencil    (&Fgrid,npoint,Even,directions,displacements),
--- a/lib/qcd/action/fermion/WilsonFermion.h
+++ b/lib/qcd/action/fermion/WilsonFermion.h
@@ -5,7 +5,7 @@ namespace Grid {

  namespace QCD {

-    class WilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
+    class WilsonFermion : public FermionOperator<LatticeFermion,LatticeGaugeField>
    {
    public:

@@ -44,7 +44,7 @@ namespace Grid {
 			int dag);

      // Constructor
-      WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass);
+      WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass);

      // DoubleStore
      void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@@ -57,7 +57,7 @@ namespace Grid {

    protected:

-      double                        mass;
+      RealD                        mass;

      GridBase                     *    _grid; 
      GridBase                     *  _cbgrid;
--- a/lib/qcd/action/fermion/FiveDimWilsonFermion.cc
+++ b/lib/qcd/action/fermion/FiveDimWilsonFermion.cc
@@ -4,18 +4,18 @@ namespace Grid {
 namespace QCD {
  
  // S-direction is INNERMOST and takes no part in the parity.
-  const std::vector<int> FiveDimWilsonFermion::directions   ({1,2,3,4, 1, 2, 3, 4});
-  const std::vector<int> FiveDimWilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
+  const std::vector<int> WilsonFermion5D::directions   ({1,2,3,4, 1, 2, 3, 4});
+  const std::vector<int> WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1});

-  int FiveDimWilsonFermion::HandOptDslash;
+  int WilsonFermion5D::HandOptDslash;

  // 5d lattice for DWF.
-  FiveDimWilsonFermion::FiveDimWilsonFermion(LatticeGaugeField &_Umu,
+  WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu,
 					   GridCartesian         &FiveDimGrid,
 					   GridRedBlackCartesian &FiveDimRedBlackGrid,
 					   GridCartesian         &FourDimGrid,
 					   GridRedBlackCartesian &FourDimRedBlackGrid,
-					   double _mass) :
+					   RealD _M5) :
  _FiveDimGrid(&FiveDimGrid),
  _FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
  _FourDimGrid(&FourDimGrid),
@@ -23,7 +23,7 @@ namespace QCD {
  Stencil    (_FiveDimGrid,npoint,Even,directions,displacements),
  StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
  StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
-  mass(_mass),
+  M5(_M5),
  Umu(_FourDimGrid),
  UmuEven(_FourDimRedBlackGrid),
  UmuOdd (_FourDimRedBlackGrid),
@@ -70,7 +70,7 @@ namespace QCD {
  pickCheckerboard(Even,UmuEven,Umu);
  pickCheckerboard(Odd ,UmuOdd,Umu);
 }
-void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
+void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
 {
  conformable(Uds._grid,GaugeGrid());
  conformable(Umu._grid,GaugeGrid());
@@ -82,60 +82,9 @@ void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const Latti
    pokeIndex<LorentzIndex>(Uds,U,mu+4);
  }
 }
-
-RealD FiveDimWilsonFermion::M(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard=in.checkerboard;
-  Dhop(in,out,DaggerNo);
-  return axpy_norm(out,5.0-M5,in,out);
-}
-RealD FiveDimWilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard=in.checkerboard;
-  Dhop(in,out,DaggerYes);
-  return axpy_norm(out,5.0-M5,in,out);
-}
-void FiveDimWilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out)
-{
-  if ( in.checkerboard == Odd ) {
-    DhopEO(in,out,DaggerNo);
-  } else {
-    DhopOE(in,out,DaggerNo);
-  }
-}
-void FiveDimWilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  if ( in.checkerboard == Odd ) {
-    DhopEO(in,out,DaggerYes);
-  } else {
-    DhopOE(in,out,DaggerYes);
-  }
-}
-void FiveDimWilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  out = (5.0-M5)*in;
-  return ;
-}
-void FiveDimWilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  Mooee(in,out);
-}
-void FiveDimWilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  out = (1.0/(5.0-M5))*in;
-  return ;
-}
-void FiveDimWilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
-{
-  out.checkerboard = in.checkerboard;
-  MooeeInv(in,out);
-}
-void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
-					LatticeDoubledGaugeField & U,
-					const LatticeFermion &in, LatticeFermion &out,int dag)
+void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
+				   LatticeDoubledGaugeField & U,
+				   const LatticeFermion &in, LatticeFermion &out,int dag)
 {
  assert((dag==DaggerNo) ||(dag==DaggerYes));

@@ -150,19 +99,21 @@ void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo
  // - 8 linear access unit stride streams per thread for Fermion for hw prefetchable.
  if ( dag == DaggerYes ) {
    if( HandOptDslash ) {
-      for(int ss=0;ss<U._grid->oSites();ss++){
-	int sU=lo.Reorder(ss);
 PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
 	for(int s=0;s<Ls;s++){
+	  //int sU=lo.Reorder(ss);
+	  int sU=ss;
 	  int sF = s+Ls*sU;
 	  DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
 	}
      }
    } else { 
-      for(int ss=0;ss<U._grid->oSites();ss++){
-	int sU=lo.Reorder(ss);
 PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
 	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
 	  int sF = s+Ls*sU;
 	  DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
 	}
@@ -170,21 +121,22 @@ PARALLEL_FOR_LOOP
    }
  } else {
    if( HandOptDslash ) {
-
 PARALLEL_FOR_LOOP
      for(int ss=0;ss<U._grid->oSites();ss++){
-	int sU=lo.Reorder(ss);
 	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
 	  int sF = s+Ls*sU;
 	  DiracOptHand::DhopSite(st,U,comm_buf,sF,sU,in,out);
 	}
      }

    } else { 
-      for(int ss=0;ss<U._grid->oSites();ss++){
-	int sU=lo.Reorder(ss);
 PARALLEL_FOR_LOOP
+      for(int ss=0;ss<U._grid->oSites();ss++){
 	for(int s=0;s<Ls;s++){
+	  //	  int sU=lo.Reorder(ss);
+	  int sU=ss;
 	  int sF = s+Ls*sU; 
 	  DiracOpt::DhopSite(st,U,comm_buf,sF,sU,in,out);
 	}
@@ -192,7 +144,7 @@ PARALLEL_FOR_LOOP
    }
  }
 }
-void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
+void WilsonFermion5D::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
 {
  conformable(in._grid,FermionRedBlackGrid());    // verifies half grid
  conformable(in._grid,out._grid); // drops the cb check
@@ -202,7 +154,7 @@ void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,

  DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
 }
-void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
+void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
 {
  conformable(in._grid,FermionRedBlackGrid());    // verifies half grid
  conformable(in._grid,out._grid); // drops the cb check
@@ -212,7 +164,7 @@ void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,

  DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
 }
-void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
+void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
 {
  conformable(in._grid,FermionGrid()); // verifies full grid
  conformable(in._grid,out._grid);
@@ -221,8 +173,14 @@ void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,in

  DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
 }
-
-}}
+void WilsonFermion5D::DW(const LatticeFermion &in, LatticeFermion &out,int dag)
+{
+  out.checkerboard=in.checkerboard;
+  Dhop(in,out,dag); // -0.5 is included
+  axpy(out,4.0-M5,in,out);
+}
+}
+}



--- a/lib/qcd/action/fermion/FiveDimWilsonFermion.h
+++ b/lib/qcd/action/fermion/FiveDimWilsonFermion.h
@@ -1,5 +1,5 @@
-#ifndef  GRID_QCD_DWF_H
-#define  GRID_QCD_DWF_H
+#ifndef  GRID_QCD_WILSON_FERMION_5D_H
+#define  GRID_QCD_WILSON_FERMION_5D_H

 namespace Grid {

@@ -14,8 +14,20 @@ namespace Grid {
    // i.e. even even contains fifth dim hopping term.
    //
    // [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ]
+    ////////////////////////////
+    //ContFrac:
+    //  Ls always odd. Rational poly deg is either Ls or Ls-1
+    //PartFrac 
+    //  Ls always odd. Rational poly deg is either Ls or Ls-1
+    //
+    //Cayley: Ls always even, Rational poly deg is Ls
+    // 
+    // Just set nrational as Ls. Forget about Ls-1 cases.
+    //
+    // Require odd Ls for cont and part frac
+    ////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////
-    class FiveDimWilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
+    class WilsonFermion5D : public FermionOperator<LatticeFermion,LatticeGaugeField>
    {
    public:
      ///////////////////////////////////////////////////////////////
@@ -26,19 +38,21 @@ namespace Grid {
      GridBase *FermionGrid(void)            { return _FiveDimGrid;}
      GridBase *FermionRedBlackGrid(void)    { return _FiveDimRedBlackGrid;}

-      // override multiply
-      virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out);
-      virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out);
+      // full checkerboard operations; leave unimplemented as abstract for now
+      //virtual RealD  M    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //virtual RealD  Mdag (const LatticeFermion &in, LatticeFermion &out)=0;

-      // half checkerboard operaions
-      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
-      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
+      // half checkerboard operations; leave unimplemented as abstract for now
+      //      virtual void   Meooe       (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MeooeDag    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   Mooee       (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeDag    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out)=0;
+      //      virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)=0;

-      // non-hermitian hopping term; half cb or both
+      // Implement hopping term non-hermitian hopping term; half cb or both
+      // Implement s-diagonal DW
+      void DW    (const LatticeFermion &in, LatticeFermion &out,int dag);
      void Dhop  (const LatticeFermion &in, LatticeFermion &out,int dag);
      void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
      void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
@@ -54,12 +68,12 @@ namespace Grid {
 			int dag);

      // Constructors
-      FiveDimWilsonFermion(LatticeGaugeField &_Umu,
+      WilsonFermion5D(LatticeGaugeField &_Umu,
 			  GridCartesian         &FiveDimGrid,
 			  GridRedBlackCartesian &FiveDimRedBlackGrid,
 			  GridCartesian         &FourDimGrid,
 			  GridRedBlackCartesian &FourDimRedBlackGrid,
-			  double _mass);
+			  double _M5);

      // DoubleStore
      void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@@ -82,7 +96,6 @@ namespace Grid {
      static const std::vector<int> displacements;

      double                        M5;
-      double                        mass;
      int Ls;

      //Defines the stencils for even and odd
--- a/lib/simd/Grid_avx.h
+++ b/lib/simd/Grid_avx.h
@@ -4,7 +4,7 @@

  Using intrinsics
 */
-// Time-stamp: <2015-05-27 12:07:15 neo>
+// Time-stamp: <2015-05-29 14:13:30 neo>
 //----------------------------------------------------------------------

 #include <immintrin.h>
@@ -261,13 +261,7 @@ namespace Optimization {
    }
    // Complex double
    inline __m256d operator()(__m256d in){
-      return _mm256_xor_pd(_mm256_addsub_pd(_mm256_setzero_pd(),in), _mm256_set1_pd(-0.f));//untested
-      /*
-	// original 
-	//      addsubps 0, inv=>0+in.v[3] 0-in.v[2], 0+in.v[1], 0-in.v[0], ...
-	__m256d tmp = _mm256_addsub_pd(_mm256_setzero_pd(),_mm256_shuffle_pd(in,in,0x5));
-	return _mm256_shuffle_pd(tmp,tmp,0x5);
-      */
+      return _mm256_xor_pd(_mm256_addsub_pd(_mm256_setzero_pd(),in), _mm256_set1_pd(-0.f));
    }
    // do not define for integer input
  };
--- a/lib/simd/Grid_vector_types.h
+++ b/lib/simd/Grid_vector_types.h
@@ -2,7 +2,7 @@
 /*! @file Grid_vector_types.h
  @brief Defines templated class Grid_simd to deal with inner vector types
 */
-// Time-stamp: <2015-05-27 12:04:06 neo>
+// Time-stamp: <2015-05-29 14:19:48 neo>
 //---------------------------------------------------------------------------
 #ifndef GRID_VECTOR_TYPES
 #define GRID_VECTOR_TYPES
@@ -55,7 +55,6 @@ namespace Grid {
  // general forms to allow for vsplat syntax
  // need explicit declaration of types when used since
  // clang cannot automatically determine the output type sometimes
-  // use decltype?
  template < class Out, class Input1, class Input2, class Operation > 
    Out binary(Input1 src_1, Input2 src_2, Operation op){
    return op(src_1, src_2);
--- a/lib/stamp-h1
+++ b/lib/stamp-h1
@@ -1 +1 @@
-timestamp for lib/Grid_config.h
+timestamp for lib/GridConfig.h
--- a/lib/stencil/Grid_lebesgue.cc
+++ b/lib/stencil/Grid_lebesgue.cc
--- a/lib/stencil/Grid_lebesgue.h
+++ b/lib/stencil/Grid_lebesgue.h
--- a/lib/stencil/Grid_stencil_common.cc
+++ b/lib/stencil/Grid_stencil_common.cc
@@ -52,8 +52,8 @@ namespace Grid {
 	// up a table containing the npoint "neighbours" and whether they 
 	// live in lattice or a comms buffer.
 	if ( !comm_dim ) {
-	  sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
-	  sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
+	  sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
+	  sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);

 	  if ( sshift[0] == sshift[1] ) {
 	    Local(point,dimension,shift,0x3);
@@ -63,8 +63,8 @@ namespace Grid {
 	  }
 	} else { // All permute extract done in comms phase prior to Stencil application
 	  //        So tables are the same whether comm_dim or splice_dim
-	  sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
-	  sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
+	  sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
+	  sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
 	  if ( sshift[0] == sshift[1] ) {
 	    Comms(point,dimension,shift,0x3);
 	  } else {
@@ -96,7 +96,7 @@ namespace Grid {
 	
 	int cb= (cbmask==0x2)? Odd : Even;
 	  
-	int sshift = _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
+	int sshift = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
 	int sx     = (x+sshift)%rd;
 	  
 	int permute_slice=0;
@@ -134,7 +134,7 @@ namespace Grid {
                                           // send to one or more remote nodes.

      int cb= (cbmask==0x2)? Odd : Even;
-      int sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
+      int sshift= _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
      
      for(int x=0;x<rd;x++){       
 	
--- a/lib/tensors/Tensor_Ta.h
+++ b/lib/tensors/Tensor_Ta.h
@@ -0,0 +1,43 @@
+#ifndef GRID_MATH_TA_H
+#define GRID_MATH_TA_H
+namespace Grid {
+
+  /////////////////////////////////////////////// 
+  // Ta function for scalar, vector, matrix
+  /////////////////////////////////////////////// 
+  inline ComplexF Ta( const ComplexF &arg){    return arg;}
+  inline ComplexD Ta( const ComplexD &arg){    return arg;}
+  inline RealF Ta( const RealF &arg){    return arg;}
+  inline RealD Ta( const RealD &arg){    return arg;}
+
+
+  template<class vtype> inline iScalar<vtype> Ta(const iScalar<vtype>&r)
+    {
+      iScalar<vtype> ret;
+      ret._internal = Ta(r._internal);
+      return ret;
+    }
+  template<class vtype,int N> inline iVector<vtype,N> Ta(const iVector<vtype,N>&r)
+    {
+      iVector<vtype,N> ret;
+      for(int i=0;i<N;i++){
+        ret._internal[i] = Ta(r._internal[i]);
+      }
+      return ret;
+    }
+  template<class vtype,int N> inline iMatrix<vtype,N> Ta(const iMatrix<vtype,N> &arg)
+    {
+      iMatrix<vtype,N> ret(arg);
+      double factor = (1/(double)N);
+      for(int c1=0;c1<N;c1++){
+	for(int c2=0;c2<N;c2++){
+	  ret._internal[c1][c2]= (ret._internal[c1][c2] - adj(arg._internal[c2][c1]));
+	  ret._internal[c1][c2] *= 0.5;
+	}}
+      //ret = (ret - adj(arg))*0.5;
+      ret -= trace(ret)*factor;
+      return ret;
+    }
+
+}
+#endif
--- a/lib/tensors/Tensor_arith.h
+++ b/lib/tensors/Tensor_arith.h
@@ -0,0 +1,11 @@
+#ifndef GRID_MATH_ARITH_H
+#define GRID_MATH_ARITH_H
+
+#include <tensors/Tensor_arith_add.h>
+#include <tensors/Tensor_arith_sub.h>
+#include <tensors/Tensor_arith_mac.h>
+#include <tensors/Tensor_arith_mul.h>
+#include <tensors/Tensor_arith_scalar.h>
+
+#endif
+
--- a/lib/tensors/Tensor_arith_add.h
+++ b/lib/tensors/Tensor_arith_add.h
@@ -45,7 +45,10 @@ namespace Grid {
  {
    for(int c2=0;c2<N;c2++){
      for(int c1=0;c1<N;c1++){
-        add(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
+	if ( c1==c2)
+	  add(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
+	else
+	  ret->_internal[c1][c2]=lhs->_internal[c1][c2];
      }}
    return;
  }
--- a/lib/tensors/Tensor_arith_mac.h
+++ b/lib/tensors/Tensor_arith_mac.h
--- a/lib/tensors/Tensor_arith_mul.h
+++ b/lib/tensors/Tensor_arith_mul.h
--- a/lib/tensors/Tensor_arith_scalar.h
+++ b/lib/tensors/Tensor_arith_scalar.h
--- a/lib/tensors/Tensor_arith_sub.h
+++ b/lib/tensors/Tensor_arith_sub.h
@@ -44,7 +44,7 @@ template<class vtype,class ltype,class rtype, int N> strong_inline void sub(iMat
                                                                     const iMatrix<rtype,N> * __restrict__ rhs){
    for(int c2=0;c2<N;c2++){
    for(int c1=0;c1<N;c1++){
-        if ( c1!=c2) {
+        if ( c1==c2) {
            sub(&ret->_internal[c1][c2],&lhs->_internal,&rhs->_internal[c1][c2]);
        } else {
            // Fails -- need unary minus. Catalogue other unops?
@@ -60,7 +60,7 @@ template<class vtype,class ltype,class rtype, int N> strong_inline void sub(iMat
                                                                     const iScalar<rtype> * __restrict__ rhs){
    for(int c2=0;c2<N;c2++){
    for(int c1=0;c1<N;c1++){
-        if ( c1!=c2)
+        if ( c1==c2)
            sub(&ret->_internal[c1][c2],&lhs->_internal[c1][c2],&rhs->_internal);
        else
            ret->_internal[c1][c2]=lhs->_internal[c1][c2];
--- a/lib/math/Grid_math_tensors.h
+++ b/lib/math/Grid_math_tensors.h
--- a/lib/tensors/Tensor_extract_merge.h
+++ b/lib/tensors/Tensor_extract_merge.h
--- a/lib/tensors/Tensor_inner.h
+++ b/lib/tensors/Tensor_inner.h
--- a/lib/tensors/Tensor_outer.h
+++ b/lib/tensors/Tensor_outer.h
--- a/lib/tensors/Tensor_peek.h
+++ b/lib/tensors/Tensor_peek.h
--- a/lib/tensors/Tensor_poke.h
+++ b/lib/tensors/Tensor_poke.h
--- a/lib/tensors/Tensor_reality.h
+++ b/lib/tensors/Tensor_reality.h
@@ -2,10 +2,6 @@
 #define GRID_MATH_REALITY_H
 namespace Grid {

-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    /////////////////////////////////////////// CONJ         ///////////////////////////////////////////
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-
 /////////////////////////////////////////////// 
 // multiply by I; make recursive.
 /////////////////////////////////////////////// 
@@ -151,6 +147,9 @@ template<class vtype,int N> inline iMatrix<vtype,N> adj(const iMatrix<vtype,N> &



+
+
+
 /////////////////////////////////////////////////////////////////
 // Can only take the real/imag part of scalar objects, since
 // lattice objects of different complex nature are non-conformable.
--- a/lib/tensors/Tensor_trace.h
+++ b/lib/tensors/Tensor_trace.h
@@ -72,5 +72,18 @@ auto traceIndex(const iMatrix<vtype,N> &arg) ->  iMatrix<decltype(traceIndex<Lev
  return ret;
 }

+// Allow to recurse if vector, but never terminate on a vector
+// trace of a different index can distribute across the vector index in a replicated way
+// but we do not trace a vector index.
+ template<int Level,class vtype,int N,typename std::enable_if< iVector<vtype, N>::TensorLevel != Level >::type * =nullptr> inline 
+auto traceIndex(const iVector<vtype,N> &arg) ->  iVector<decltype(traceIndex<Level>(arg._internal[0])),N> 
+{
+  iVector<decltype(traceIndex<Level>(arg._internal[0])),N> ret;
+  for(int i=0;i<N;i++){
+    ret._internal[i] = traceIndex<Level>(arg._internal[i]);
+  }
+  return ret;
+}
+
 }
 #endif
--- a/Show More
+++ b/Show More