PartialFraction Hw with Zolo and Tanh approx converged under CG and passed EO breakdown

and hermiticity tests.

INSTALL

@@ -1 +1 @@
-/usr/share/automake-1.14/INSTALL
+/opt/local/share/automake-1.15/INSTALL

configure

@@ -2574,7 +2574,7 @@ test -n "$target_alias" &&
     NONENONEs,x,x, &&
   program_prefix=${target_alias}-
 
-am__api_version='1.14'
+am__api_version='1.15'
 
 # Find a good install program.  We prefer a C program (faster),
 # so one script is as good as another.  But avoid the broken or
@@ -2746,8 +2746,8 @@ test "$program_suffix" != NONE &&
 ac_script='s/[\\$]/&&/g;s/;s,x,x,$//'
 program_transform_name=`$as_echo "$program_transform_name" | sed "$ac_script"`
 
-# expand $ac_aux_dir to an absolute path
+# Expand $ac_aux_dir to an absolute path.
-am_aux_dir=`cd $ac_aux_dir && pwd`
+am_aux_dir=`cd "$ac_aux_dir" && pwd`
 
 if test x"${MISSING+set}" != xset; then
   case $am_aux_dir in
@@ -2766,7 +2766,7 @@ else
 $as_echo "$as_me: WARNING: 'missing' script is too old or missing" >&2;}
 fi
 
-if test x"${install_sh}" != xset; then
+if test x"${install_sh+set}" != xset; then
   case $am_aux_dir in
   *\ * | *\	*)
     install_sh="\${SHELL} '$am_aux_dir/install-sh'" ;;
@@ -3094,8 +3094,8 @@ MAKEINFO=${MAKEINFO-"${am_missing_run}makeinfo"}
 # <http://lists.gnu.org/archive/html/automake/2012-07/msg00014.html>
 mkdir_p='$(MKDIR_P)'
 
-# We need awk for the "check" target.  The system "awk" is bad on
+# We need awk for the "check" target (and possibly the TAP driver).  The
-# some platforms.
+# system "awk" is bad on some platforms.
 # Always define AMTAR for backward compatibility.  Yes, it's still used
 # in the wild :-(  We should find a proper way to deprecate it ...
 AMTAR='$${TAR-tar}'
@@ -3154,6 +3154,7 @@ END
 fi
 
 
+
 ac_config_headers="$ac_config_headers lib/GridConfig.h"
 
 # Check whether --enable-silent-rules was given.
@@ -6656,9 +6657,6 @@ fi
 
 
 
-
-
-
 # Check whether --enable-simd was given.
 if test "${enable_simd+set}" = set; then :
   enableval=$enable_simd; \
@@ -6673,21 +6671,21 @@ supported=no
 case ${ac_SIMD} in
      SSE4)
        echo Configuring for SSE4
-       if test x"$ax_cv_support_ssse3_ext" = x"yes"; then
+
 $as_echo "#define SSE4 1" >>confdefs.h
 
-       supported=yes
+       if test x"$ax_cv_support_ssse3_ext" = x"yes"; then           supported=yes
        else
-	{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: Your processor does not support SSE4 instructions" >&5
+  	{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: Your processor does not support SSE4 instructions" >&5
 $as_echo "$as_me: WARNING: Your processor does not support SSE4 instructions" >&2;}
        fi
      ;;
      AVX)
        echo Configuring for AVX
-       if test x"$ax_cv_support_avx_ext" = x"yes"; then
+
 $as_echo "#define AVX1 1" >>confdefs.h
 
-       supported=yes
+       if test x"$ax_cv_support_avx_ext" = x"yes"; then         supported=yes
        else
        	{ $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: Your processor does not support AVX instructions" >&5
 $as_echo "$as_me: WARNING: Your processor does not support AVX instructions" >&2;}
@@ -6695,10 +6693,10 @@ $as_echo "$as_me: WARNING: Your processor does not support AVX instructions" >&2
      ;;
      AVX2)
        echo Configuring for AVX2
-       if test x"$ax_cv_support_avx2_ext" = x"yes"; then
+
 $as_echo "#define AVX2 1" >>confdefs.h
 
-       supported=yes
+       if test x"$ax_cv_support_avx2_ext" = x"yes"; then         supported=yes
        else
        { $as_echo "$as_me:${as_lineno-$LINENO}: WARNING: Your processor does not support AVX2 instructions" >&5
 $as_echo "$as_me: WARNING: Your processor does not support AVX2 instructions" >&2;}

configure.ac

@@ -66,9 +66,6 @@ Please install or provide the correct path to your installation
 Info at: http://www.mpfr.org/)])
 
 
-
-
-
 AC_ARG_ENABLE([simd],[AC_HELP_STRING([--enable-simd=SSE4|AVX|AVX2|AVX512|MIC],\
 	[Select instructions to be SSE4.0, AVX 1.0, AVX 2.0+FMA, AVX 512, MIC])],\
 	[ac_SIMD=${enable_simd}],[ac_SIMD=AVX2])
@@ -78,17 +75,17 @@ supported=no
 case ${ac_SIMD} in
      SSE4)
        echo Configuring for SSE4
-       if test x"$ax_cv_support_ssse3_ext" = x"yes"; then  dnl minimal support for SSE4
        AC_DEFINE([SSE4],[1],[SSE4] )
-       supported=yes
+       if test x"$ax_cv_support_ssse3_ext" = x"yes"; then  dnl minimal support for SSE4
+         supported=yes
        else
-	AC_MSG_WARN([Your processor does not support SSE4 instructions])
+  	AC_MSG_WARN([Your processor does not support SSE4 instructions])
        fi
      ;;
      AVX)
        echo Configuring for AVX
-       if test x"$ax_cv_support_avx_ext" = x"yes"; then  dnl minimal support for AVX
        AC_DEFINE([AVX1],[1],[AVX] )
+       if test x"$ax_cv_support_avx_ext" = x"yes"; then  dnl minimal support for AVX
        supported=yes			  
        else
        	AC_MSG_WARN([Your processor does not support AVX instructions])
@@ -96,8 +93,8 @@ case ${ac_SIMD} in
      ;;
      AVX2)
        echo Configuring for AVX2
-       if test x"$ax_cv_support_avx2_ext" = x"yes"; then  dnl minimal support for AVX2
        AC_DEFINE([AVX2],[1],[AVX2] )
+       if test x"$ax_cv_support_avx2_ext" = x"yes"; then  dnl minimal support for AVX2
        supported=yes
        else
        AC_MSG_WARN([Your processor does not support AVX2 instructions])

lib/Make.inc

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

lib/algorithms/approx/Zolotarev.cc

@@ -293,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* grid_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;
@@ -414,7 +414,19 @@ zolotarev_data* grid_zolotarev(PRECISION epsilon, int n, int type) {
   return zd;
 }
 
-zolotarev_data* grid_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;

lib/algorithms/approx/Zolotarev.h

@@ -77,8 +77,9 @@ typedef struct {
  * zolotarev_data structure. The arguments must satisfy the constraints that
  * epsilon > 0, n > 0, and type = 0 or 1. */
 
-ZOLOTAREV_DATA* grid_higham(PRECISION epsilon, int n) ;
+ZOLOTAREV_DATA* higham(PRECISION epsilon, int n) ;
-ZOLOTAREV_DATA* grid_zolotarev(PRECISION epsilon, int n, int type);
+ZOLOTAREV_DATA* zolotarev(PRECISION epsilon, int n, int type);
+void zolotarev_free(zolotarev_data *zdata);
 #endif
 
 #ifdef __cplusplus

lib/qcd/action/Actions.h

@@ -62,6 +62,8 @@
 // Partial fraction
 //////////////////////
 #include <qcd/action/fermion/PartialFractionFermion5D.h>
+#include <qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h>
+#include <qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h>
 
 
     // Chroma interface defining FermionAction

lib/qcd/action/fermion/ContinuedFractionFermion5D.cc

@@ -9,6 +9,16 @@ namespace Grid {
     }
     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);
@@ -29,7 +39,7 @@ namespace Grid {
 
 
       ZoloHiInv =1.0/zolo_hi;
-      double dw_diag = (4.0-M5)*ZoloHiInv;
+      dw_diag = (4.0-M5)*ZoloHiInv;
     
       See.resize(Ls);
       Aee.resize(Ls);
@@ -105,8 +115,6 @@ namespace Grid {
     }
     void   ContinuedFractionFermion5D::Mooee       (const LatticeFermion &psi, LatticeFermion &chi)
     {
-      double dw_diag = (4.0-M5)*ZoloHiInv;
-    
       int sign=1;
       for(int s=0;s<Ls;s++){
 	if ( s==0 ) {

lib/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -37,9 +37,8 @@ namespace Grid {
       void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
       void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);;
 
-      Approx::zolotarev_data *zdata;
-
       // Cont frac
+      RealD dw_diag;
       RealD mass;
       RealD R;
       RealD ZoloHiInv;

lib/qcd/action/fermion/DomainWallFermion.h

@@ -29,13 +29,14 @@ namespace Grid {
       {
 	RealD eps = 1.0;
 
-	Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
+	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);
       }
 
     };

lib/qcd/action/fermion/MobiusFermion.h

@@ -31,11 +31,13 @@ namespace Grid {
 	RealD eps = 1.0;
 
 	std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
-	Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
+	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);
  
       }
 

lib/qcd/action/fermion/MobiusZolotarevFermion.h

@@ -31,7 +31,7 @@ namespace Grid {
       {
 	RealD eps = lo/hi;
 
-	Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,this->Ls,0);// eps is ignored for higham
+	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;
@@ -39,6 +39,7 @@ namespace Grid {
 	// Call base setter
 	this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);
  
+	Approx::zolotarev_free(zdata);
       }
 
     };

lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h

@@ -30,8 +30,9 @@ namespace Grid {
 	{
 	  assert((Ls&0x1)==1); // Odd Ls required
 	  int nrational=Ls-1;// Even rational order
-	  zdata = Approx::grid_higham(1.0,nrational);// eps is ignored for higham
+	  Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham
 	  SetCoefficientsTanh(zdata,scale);
+	  Approx::zolotarev_free(zdata);
 	}
     };
   }

lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h

@@ -30,12 +30,12 @@ namespace Grid {
 	{
 	  assert((Ls&0x1)==1); // Odd Ls required
 
-	  int nrational=Ls-1;// Even rational order
+	  int nrational=Ls;// Odd rational order
 	  RealD eps = lo/hi;
 
-	  Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,nrational,0);
+	  Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0);
-
 	  SetCoefficientsZolotarev(hi,zdata);
+	  Approx::zolotarev_free(zdata);
 
 	}
     };

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

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

lib/qcd/action/fermion/PartialFractionFermion5D.cc

@@ -1 +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);
+
+    }
+ 
+ }
+}
 

lib/qcd/action/fermion/PartialFractionFermion5D.h

@@ -9,6 +9,13 @@ namespace Grid {
     {
     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);
@@ -21,16 +28,7 @@ namespace Grid {
       virtual void   MooeeInv    (const LatticeFermion &in, LatticeFermion &out);
       virtual void   MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
 
-    private:
+      virtual void   Instantiatable(void) =0; // ensure no make-eee
-
-      virtual void PartialFractionCoefficients(void);
-
-      Approx::zolotarev_data *zdata;
-
-      // Part frac
-      double R;
-      std::vector<double> p; 
-      std::vector<double> q;
 
       // Constructors
       PartialFractionFermion5D(LatticeGaugeField &_Umu,
@@ -40,6 +38,20 @@ namespace Grid {
 				    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;
+
     };
 
 

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -1,5 +1,5 @@
-#ifndef  GRID_QCD_DWF_H
+#ifndef  GRID_QCD_WILSON_FERMION_5D_H
-#define  GRID_QCD_DWF_H
+#define  GRID_QCD_WILSON_FERMION_5D_H
 
 namespace Grid {
 

tests/Test_contfrac_cg.cc

@@ -82,6 +82,16 @@ int main (int argc, char ** argv)
   OverlapWilsonContFracZolotarevFermion Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
   TestCGinversions<OverlapWilsonContFracZolotarevFermion>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+
+  std::cout <<"OverlapWilsonPartialFractionTanhFermion  test"<<std::endl;
+  OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+  TestCGinversions<OverlapWilsonPartialFractionTanhFermion>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
+  std::cout <<"OverlapWilsonPartialFractionZolotarevFermion  test"<<std::endl;
+  OverlapWilsonPartialFractionZolotarevFermion Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
+  TestCGinversions<OverlapWilsonPartialFractionZolotarevFermion>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
+
   Grid_finalize();
 }
 template<class What> 

tests/Test_contfrac_even_odd.cc

@@ -57,6 +57,14 @@ int main (int argc, char ** argv)
   OverlapWilsonContFracZolotarevFermion Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
   TestWhat<OverlapWilsonContFracZolotarevFermion>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout <<"OverlapWilsonPartialFractionTanhFermion  test"<<std::endl;
+  OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+  TestWhat<OverlapWilsonPartialFractionTanhFermion>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
+  std::cout <<"OverlapWilsonPartialFractionZolotarevFermion  test"<<std::endl;
+  OverlapWilsonPartialFractionZolotarevFermion Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
+  TestWhat<OverlapWilsonPartialFractionZolotarevFermion>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+
   Grid_finalize();
 }