Merge branch 'develop' into feature/gpt

Grid/Grid_Eigen_Dense.h

@@ -22,8 +22,18 @@
 #undef __CUDACC__
 #undef __CUDA_ARCH__
 #define __NVCC__REDEFINE__
+#endif 
+
+/* SYCL save and restore compile environment*/
+#ifdef __SYCL_DEVICE_ONLY__  
+#pragma push
+#pragma push_macro("__SYCL_DEVICE_ONLY__")
+#undef __SYCL_DEVICE_ONLY__
+#undef EIGEN_USE_SYCL
+#define EIGEN_DONT_VECTORIZE
 #endif
 
+
 #include <Grid/Eigen/Dense>
 #include <Grid/Eigen/unsupported/CXX11/Tensor>
 
@@ -35,7 +45,14 @@
 #pragma pop
 #endif
 
+/*SYCL restore*/
+#ifdef __SYCL__REDEFINE__
+#pragma pop_macro("__SYCL_DEVICE_ONLY__")
+#pragma pop
+#endif
+
 #if defined __GNUC__
 #pragma GCC diagnostic pop
 #endif
 
+

Grid/algorithms/Algorithms.h

@@ -39,14 +39,18 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/approx/Remez.h>
 #include <Grid/algorithms/approx/MultiShiftFunction.h>
 #include <Grid/algorithms/approx/Forecast.h>
+#include <Grid/algorithms/approx/RemezGeneral.h>
+#include <Grid/algorithms/approx/ZMobius.h>
 
 #include <Grid/algorithms/iterative/Deflation.h>
 #include <Grid/algorithms/iterative/ConjugateGradient.h>
+#include <Grid/algorithms/iterative/BiCGSTAB.h>
 #include <Grid/algorithms/iterative/ConjugateResidual.h>
 #include <Grid/algorithms/iterative/NormalEquations.h>
 #include <Grid/algorithms/iterative/SchurRedBlack.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMultiShift.h>
 #include <Grid/algorithms/iterative/ConjugateGradientMixedPrec.h>
+#include <Grid/algorithms/iterative/BiCGSTABMixedPrec.h>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
 #include <Grid/algorithms/iterative/MinimalResidual.h>

Grid/algorithms/LinearOperator.h

@@ -334,6 +334,132 @@ public:
 	return axpy_norm(out,-1.0,tmp,in);
       }
     };
+
+    template<class Field>
+    class NonHermitianSchurOperatorBase :  public LinearOperatorBase<Field> 
+    {
+      public:
+        virtual RealD Mpc      (const Field& in, Field& out) = 0;
+        virtual RealD MpcDag   (const Field& in, Field& out) = 0;
+        virtual void  MpcDagMpc(const Field& in, Field& out, RealD& ni, RealD& no) {
+          Field tmp(in.Grid());
+          tmp.Checkerboard() = in.Checkerboard();
+	        ni = Mpc(in,tmp);
+	        no = MpcDag(tmp,out);
+        }
+        virtual void HermOpAndNorm(const Field& in, Field& out, RealD& n1, RealD& n2) {
+          assert(0);
+        }
+        virtual void HermOp(const Field& in, Field& out) {
+          assert(0);
+        }
+        void Op(const Field& in, Field& out) {
+          Mpc(in, out);
+        }
+        void AdjOp(const Field& in, Field& out) { 
+          MpcDag(in, out);
+        }
+        // Support for coarsening to a multigrid
+        void OpDiag(const Field& in, Field& out) {
+          assert(0); // must coarsen the unpreconditioned system
+        }
+        void OpDir(const Field& in, Field& out, int dir, int disp) {
+          assert(0);
+        }
+    };
+
+    template<class Matrix, class Field>
+    class NonHermitianSchurDiagMooeeOperator :  public NonHermitianSchurOperatorBase<Field> 
+    {
+      public:
+        Matrix& _Mat;
+        NonHermitianSchurDiagMooeeOperator(Matrix& Mat): _Mat(Mat){};
+        virtual RealD Mpc(const Field& in, Field& out) {
+          Field tmp(in.Grid());
+          tmp.Checkerboard() = !in.Checkerboard();
+
+  	      _Mat.Meooe(in, tmp);
+	        _Mat.MooeeInv(tmp, out);
+	        _Mat.Meooe(out, tmp);
+
+	        _Mat.Mooee(in, out);
+	
+          return axpy_norm(out, -1.0, tmp, out);
+        }
+        virtual RealD MpcDag(const Field& in, Field& out) {
+	        Field tmp(in.Grid());
+
+	        _Mat.MeooeDag(in, tmp);
+          _Mat.MooeeInvDag(tmp, out);
+	        _Mat.MeooeDag(out, tmp);
+
+	        _Mat.MooeeDag(in, out);
+	
+          return axpy_norm(out, -1.0, tmp, out);
+      }
+    };
+    
+    template<class Matrix,class Field>
+    class NonHermitianSchurDiagOneOperator : public NonHermitianSchurOperatorBase<Field> 
+    {
+      protected:
+        Matrix &_Mat;
+    
+      public:
+        NonHermitianSchurDiagOneOperator (Matrix& Mat): _Mat(Mat){};
+        virtual RealD Mpc(const Field& in, Field& out) {
+	        Field tmp(in.Grid());
+
+	        _Mat.Meooe(in, out);
+	        _Mat.MooeeInv(out, tmp);
+	        _Mat.Meooe(tmp, out);
+	        _Mat.MooeeInv(out, tmp);
+
+	        return axpy_norm(out, -1.0, tmp, in);
+        }
+        virtual RealD MpcDag(const Field& in, Field& out) {
+	        Field tmp(in.Grid());
+
+	        _Mat.MooeeInvDag(in, out);
+	        _Mat.MeooeDag(out, tmp);
+	        _Mat.MooeeInvDag(tmp, out);
+	        _Mat.MeooeDag(out, tmp);
+
+	        return axpy_norm(out, -1.0, tmp, in);
+        }
+    };
+
+    template<class Matrix, class Field>
+    class NonHermitianSchurDiagTwoOperator : public NonHermitianSchurOperatorBase<Field> 
+    {
+      protected:
+        Matrix& _Mat;
+    
+      public:
+        NonHermitianSchurDiagTwoOperator(Matrix& Mat): _Mat(Mat){};
+
+        virtual RealD Mpc(const Field& in, Field& out) {
+          Field tmp(in.Grid());
+
+	        _Mat.MooeeInv(in, out);
+	        _Mat.Meooe(out, tmp);
+	        _Mat.MooeeInv(tmp, out);
+	        _Mat.Meooe(out, tmp);
+
+	        return axpy_norm(out, -1.0, tmp, in);
+        }
+        virtual RealD MpcDag(const Field& in, Field& out) {
+	        Field tmp(in.Grid());
+
+          _Mat.MeooeDag(in, out);
+          _Mat.MooeeInvDag(out, tmp);
+          _Mat.MeooeDag(tmp, out);
+          _Mat.MooeeInvDag(out, tmp);
+
+          return axpy_norm(out, -1.0, tmp, in);
+        }
+    };
+
     ///////////////////////////////////////////////////////////////////////////////////////////////////
     // Left  handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) psi = eta  -->  ( 1 - Moo^-1 Moe Mee^-1 Meo ) psi = Moo^-1 eta
     // Right handed Moo^-1 ; (Moo - Moe Mee^-1 Meo) Moo^-1 Moo psi = eta  -->  ( 1 - Moe Mee^-1 Meo Moo^-1) phi=eta ; psi = Moo^-1 phi

Grid/algorithms/approx/RemezGeneral.cc

@@ -0,0 +1,473 @@
+#include<math.h>
+#include<stdio.h>
+#include<stdlib.h>
+#include<string>
+#include<iostream>
+#include<iomanip>
+#include<cassert>
+
+#include<Grid/algorithms/approx/RemezGeneral.h>
+
+
+// Constructor
+AlgRemezGeneral::AlgRemezGeneral(double lower, double upper, long precision,
+				 bigfloat (*f)(bigfloat x, void *data), void *data): f(f), 
+										     data(data), 
+										     prec(precision),
+										     apstrt(lower), apend(upper), apwidt(upper - lower),
+										     n(0), d(0), pow_n(0), pow_d(0)
+{
+  bigfloat::setDefaultPrecision(prec);
+
+  std::cout<<"Approximation bounds are ["<<apstrt<<","<<apend<<"]\n";
+  std::cout<<"Precision of arithmetic is "<<precision<<std::endl;
+}
+
+//Determine the properties of the numerator and denominator polynomials
+void AlgRemezGeneral::setupPolyProperties(int num_degree, int den_degree, PolyType num_type_in, PolyType den_type_in){
+  pow_n = num_degree;
+  pow_d = den_degree;
+
+  if(pow_n % 2 == 0 && num_type_in == PolyType::Odd) assert(0);
+  if(pow_n % 2 == 1 && num_type_in == PolyType::Even) assert(0);
+
+  if(pow_d % 2 == 0 && den_type_in == PolyType::Odd) assert(0);
+  if(pow_d % 2 == 1 && den_type_in == PolyType::Even) assert(0);
+
+  num_type = num_type_in;
+  den_type = den_type_in;
+
+  num_pows.resize(pow_n+1);
+  den_pows.resize(pow_d+1);
+
+  int n_in = 0;
+  bool odd = num_type == PolyType::Full || num_type == PolyType::Odd;
+  bool even = num_type == PolyType::Full || num_type == PolyType::Even;
+  for(int i=0;i<=pow_n;i++){
+    num_pows[i] = -1;
+    if(i % 2 == 0 && even) num_pows[i] = n_in++;
+    if(i % 2 == 1 && odd) num_pows[i] = n_in++;
+  }
+
+  std::cout << n_in << " terms in numerator" << std::endl;
+  --n_in; //power is 1 less than the number of terms, eg  pow=1   a x^1  + b x^0
+
+  int d_in = 0;
+  odd = den_type == PolyType::Full || den_type == PolyType::Odd;
+  even = den_type == PolyType::Full || den_type == PolyType::Even;
+  for(int i=0;i<=pow_d;i++){
+    den_pows[i] = -1;
+    if(i % 2 == 0 && even) den_pows[i] = d_in++;
+    if(i % 2 == 1 && odd) den_pows[i] = d_in++;
+  }
+
+  std::cout << d_in << " terms in denominator" << std::endl;
+  --d_in;
+
+  n = n_in;
+  d = d_in;
+}
+
+//Setup algorithm
+void AlgRemezGeneral::reinitializeAlgorithm(){
+  spread = 1.0e37;
+  iter = 0;
+
+  neq = n + d + 1; //not +2 because highest-power term in denominator is fixed to 1
+
+  param.resize(neq);
+  yy.resize(neq+1);
+
+  //Initialize linear equation temporaries
+  A.resize(neq*neq);
+  B.resize(neq);
+  IPS.resize(neq);
+
+  //Initialize maximum and minimum errors
+  xx.resize(neq+2);
+  mm.resize(neq+1);
+  initialGuess();
+
+  //Initialize search steps
+  step.resize(neq+1);
+  stpini();
+}
+
+double AlgRemezGeneral::generateApprox(const int num_degree, const int den_degree, 
+				       const PolyType num_type_in, const PolyType den_type_in, 
+				       const double _tolerance, const int report_freq){
+  //Setup the properties of the polynomial
+  setupPolyProperties(num_degree, den_degree, num_type_in, den_type_in);
+
+  //Setup the algorithm
+  reinitializeAlgorithm();
+
+  bigfloat tolerance = _tolerance;
+
+  //Iterate until convergance
+  while (spread > tolerance) { 
+    if (iter++ % report_freq==0)
+      std::cout<<"Iteration " <<iter-1<<" spread "<<(double)spread<<" delta "<<(double)delta << std::endl; 
+
+    equations();
+    if (delta < tolerance) {
+      std::cout<<"Iteration " << iter-1 << " delta too small (" << delta << "<" << tolerance << "), try increasing precision\n";
+      assert(0);
+    };    
+    assert( delta>= tolerance );
+
+    search();
+  }
+
+  int sign;
+  double error = (double)getErr(mm[0],&sign);
+  std::cout<<"Converged at "<<iter<<" iterations; error = "<<error<<std::endl;
+
+  // Return the maximum error in the approximation
+  return error;
+}
+
+
+// Initial values of maximal and minimal errors
+void AlgRemezGeneral::initialGuess(){
+  // Supply initial guesses for solution points
+  long ncheb = neq;			// Degree of Chebyshev error estimate
+
+  // Find ncheb+1 extrema of Chebyshev polynomial
+  bigfloat a = ncheb;
+  bigfloat r;
+
+  mm[0] = apstrt;
+  for (long i = 1; i < ncheb; i++) {
+    r = 0.5 * (1 - cos((M_PI * i)/(double) a));
+    //r *= sqrt_bf(r);
+    r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
+    mm[i] = apstrt + r * apwidt;
+  }
+  mm[ncheb] = apend;
+
+  a = 2.0 * ncheb;
+  for (long i = 0; i <= ncheb; i++) {
+    r = 0.5 * (1 - cos(M_PI * (2*i+1)/(double) a));
+    //r *= sqrt_bf(r); // Squeeze to low end of interval
+    r = (exp((double)r)-1.0)/(exp(1.0)-1.0);
+    xx[i] = apstrt + r * apwidt;
+  }
+}
+
+// Initialise step sizes
+void AlgRemezGeneral::stpini(){
+  xx[neq+1] = apend;
+  delta = 0.25;
+  step[0] = xx[0] - apstrt;
+  for (int i = 1; i < neq; i++) step[i] = xx[i] - xx[i-1];
+  step[neq] = step[neq-1];
+}
+
+// Search for error maxima and minima
+void AlgRemezGeneral::search(){
+  bigfloat a, q, xm, ym, xn, yn, xx1;
+  int emsign, ensign, steps;
+
+  int meq = neq + 1;
+
+  bigfloat eclose = 1.0e30;
+  bigfloat farther = 0l;
+
+  bigfloat xx0 = apstrt;
+
+  for (int i = 0; i < meq; i++) {
+    steps = 0;
+    xx1 = xx[i]; // Next zero
+    if (i == meq-1) xx1 = apend;
+    xm = mm[i];
+    ym = getErr(xm,&emsign);
+    q = step[i];
+    xn = xm + q;
+    if (xn < xx0 || xn >= xx1) {	// Cannot skip over adjacent boundaries
+      q = -q;
+      xn = xm;
+      yn = ym;
+      ensign = emsign;
+    } else {
+      yn = getErr(xn,&ensign);
+      if (yn < ym) {
+	q = -q;
+	xn = xm;
+	yn = ym;
+	ensign = emsign;
+      }
+    }
+  
+    while(yn >= ym) {		// March until error becomes smaller.
+      if (++steps > 10)
+      	break;
+      
+      ym = yn;
+      xm = xn;
+      emsign = ensign;
+      a = xm + q;
+      if (a == xm || a <= xx0 || a >= xx1)
+	break;// Must not skip over the zeros either side.      
+
+      xn = a;
+      yn = getErr(xn,&ensign);
+    }
+
+    mm[i] = xm;			// Position of maximum
+    yy[i] = ym;			// Value of maximum
+
+    if (eclose > ym) eclose = ym;
+    if (farther < ym) farther = ym;
+
+    xx0 = xx1; // Walk to next zero.
+  } // end of search loop
+
+  q = (farther - eclose);	// Decrease step size if error spread increased
+
+  if (eclose != 0.0) q /= eclose; // Relative error spread
+
+  if (q >= spread)
+    delta *= 0.5; // Spread is increasing; decrease step size
+  
+  spread = q;
+
+  for (int i = 0; i < neq; i++) {
+    q = yy[i+1];
+    if (q != 0.0) q = yy[i] / q  - (bigfloat)1l;
+    else q = 0.0625;
+    if (q > (bigfloat)0.25) q = 0.25;
+    q *= mm[i+1] - mm[i];
+    step[i] = q * delta;
+  }
+  step[neq] = step[neq-1];
+  
+  for (int i = 0; i < neq; i++) {	// Insert new locations for the zeros.
+    xm = xx[i] - step[i];
+
+    if (xm <= apstrt)
+      continue;
+
+    if (xm >= apend)
+      continue;
+
+    if (xm <= mm[i])
+      xm = (bigfloat)0.5 * (mm[i] + xx[i]);    
+
+    if (xm >= mm[i+1])
+      xm = (bigfloat)0.5 * (mm[i+1] + xx[i]);
+    
+    xx[i] = xm;
+  }
+}
+
+// Solve the equations
+void AlgRemezGeneral::equations(){
+  bigfloat x, y, z;
+  bigfloat *aa;
+  
+  for (int i = 0; i < neq; i++) {	// set up the equations for solution by simq()
+    int ip = neq * i;		// offset to 1st element of this row of matrix
+    x = xx[i];			// the guess for this row
+    y = func(x);		// right-hand-side vector
+
+    z = (bigfloat)1l;
+    aa = A.data()+ip;
+    int t = 0;
+    for (int j = 0; j <= pow_n; j++) {
+      if(num_pows[j] != -1){ *aa++ = z; t++; }
+      z *= x;
+    }
+    assert(t == n+1);
+
+    z = (bigfloat)1l;
+    t = 0;
+    for (int j = 0; j < pow_d; j++) {
+      if(den_pows[j] != -1){ *aa++ = -y * z; t++; }
+      z *= x;
+    }
+    assert(t == d);
+
+    B[i] = y * z;		// Right hand side vector
+  }
+
+  // Solve the simultaneous linear equations.
+  if (simq()){
+    std::cout<<"simq failed\n";
+    exit(0);
+  }
+}
+
+
+// Evaluate the rational form P(x)/Q(x) using coefficients
+// from the solution vector param
+bigfloat AlgRemezGeneral::approx(const bigfloat x) const{
+  // Work backwards toward the constant term.
+  int c = n;
+  bigfloat yn = param[c--];		// Highest order numerator coefficient
+  for (int i = pow_n-1; i >= 0; i--) yn = x * yn  +  (num_pows[i] != -1 ? param[c--] : bigfloat(0l));  
+
+  c = n+d;
+  bigfloat yd = 1l; //Highest degree coefficient is 1.0
+  for (int i = pow_d-1; i >= 0; i--) yd = x * yd  +  (den_pows[i] != -1 ? param[c--] : bigfloat(0l)); 
+
+  return(yn/yd);
+}
+
+// Compute size and sign of the approximation error at x
+bigfloat AlgRemezGeneral::getErr(bigfloat x, int *sign) const{
+  bigfloat f = func(x);
+  bigfloat e = approx(x) - f;
+  if (f != 0) e /= f;
+  if (e < (bigfloat)0.0) {
+    *sign = -1;
+    e = -e;
+  }
+  else *sign = 1;
+  
+  return(e);
+}
+
+// Solve the system AX=B
+int AlgRemezGeneral::simq(){
+
+  int ip, ipj, ipk, ipn;
+  int idxpiv;
+  int kp, kp1, kpk, kpn;
+  int nip, nkp;
+  bigfloat em, q, rownrm, big, size, pivot, sum;
+  bigfloat *aa;
+  bigfloat *X = param.data();
+
+  int n = neq;
+  int nm1 = n - 1;
+  // Initialize IPS and X
+  
+  int ij = 0;
+  for (int i = 0; i < n; i++) {
+    IPS[i] = i;
+    rownrm = 0.0;
+    for(int j = 0; j < n; j++) {
+      q = abs_bf(A[ij]);
+      if(rownrm < q) rownrm = q;
+      ++ij;
+    }
+    if (rownrm == (bigfloat)0l) {
+      std::cout<<"simq rownrm=0\n";
+      return(1);
+    }
+    X[i] = (bigfloat)1.0 / rownrm;
+  }
+  
+  for (int k = 0; k < nm1; k++) {
+    big = 0.0;
+    idxpiv = 0;
+    
+    for (int i = k; i < n; i++) {
+      ip = IPS[i];
+      ipk = n*ip + k;
+      size = abs_bf(A[ipk]) * X[ip];
+      if (size > big) {
+	big = size;
+	idxpiv = i;
+      }
+    }
+    
+    if (big == (bigfloat)0l) {
+      std::cout<<"simq big=0\n";
+      return(2);
+    }
+    if (idxpiv != k) {
+      int j = IPS[k];
+      IPS[k] = IPS[idxpiv];
+      IPS[idxpiv] = j;
+    }
+    kp = IPS[k];
+    kpk = n*kp + k;
+    pivot = A[kpk];
+    kp1 = k+1;
+    for (int i = kp1; i < n; i++) {
+      ip = IPS[i];
+      ipk = n*ip + k;
+      em = -A[ipk] / pivot;
+      A[ipk] = -em;
+      nip = n*ip;
+      nkp = n*kp;
+      aa = A.data()+nkp+kp1;
+      for (int j = kp1; j < n; j++) {
+	ipj = nip + j;
+	A[ipj] = A[ipj] + em * *aa++;
+      }
+    }
+  }
+  kpn = n * IPS[n-1] + n - 1;	// last element of IPS[n] th row
+  if (A[kpn] == (bigfloat)0l) {
+    std::cout<<"simq A[kpn]=0\n";
+    return(3);
+  }
+
+  
+  ip = IPS[0];
+  X[0] = B[ip];
+  for (int i = 1; i < n; i++) {
+    ip = IPS[i];
+    ipj = n * ip;
+    sum = 0.0;
+    for (int j = 0; j < i; j++) {
+      sum += A[ipj] * X[j];
+      ++ipj;
+    }
+    X[i] = B[ip] - sum;
+  }
+  
+  ipn = n * IPS[n-1] + n - 1;
+  X[n-1] = X[n-1] / A[ipn];
+  
+  for (int iback = 1; iback < n; iback++) {
+    //i goes (n-1),...,1
+    int i = nm1 - iback;
+    ip = IPS[i];
+    nip = n*ip;
+    sum = 0.0;
+    aa = A.data()+nip+i+1;
+    for (int j= i + 1; j < n; j++) 
+      sum += *aa++ * X[j];
+    X[i] = (X[i] - sum) / A[nip+i];
+  }
+  
+  return(0);
+}
+
+void AlgRemezGeneral::csv(std::ostream & os) const{
+  os << "Numerator" << std::endl;
+  for(int i=0;i<=pow_n;i++){
+    os << getCoeffNum(i) << "*x^" << i;
+    if(i!=pow_n) os << " + ";
+  }
+  os << std::endl;
+
+  os << "Denominator" << std::endl;
+  for(int i=0;i<=pow_d;i++){
+    os << getCoeffDen(i) << "*x^" << i;
+    if(i!=pow_d) os << " + ";
+  }
+  os << std::endl;
+
+  //For a true minimax solution the errors should all be equal and the signs should oscillate +-+-+- etc
+  int sign;
+  os << "Errors at maxima: coordinate, error, (sign)" << std::endl;
+  for(int i=0;i<neq+1;i++){ 
+    os << mm[i] << " " << getErr(mm[i],&sign) << " (" << sign << ")" << std::endl;
+  }
+
+  os << "Scan over range:" << std::endl;
+  int npt = 60;
+  bigfloat dlt = (apend - apstrt)/bigfloat(npt-1);
+
+  for (bigfloat x=apstrt; x<=apend; x = x + dlt) {
+    double f = evaluateFunc(x);
+    double r = evaluateApprox(x);
+    os<< x<<","<<r<<","<<f<<","<<r-f<<std::endl;
+  }
+  return;
+}

Grid/algorithms/approx/RemezGeneral.h

@@ -0,0 +1,170 @@
+/*
+  C.Kelly Jan 2020 based on implementation by M. Clark May 2005
+
+  AlgRemezGeneral is an implementation of the Remez algorithm for approximating an arbitrary function by a rational polynomial 
+  It includes optional restriction to odd/even polynomials for the numerator and/or denominator
+*/
+
+#ifndef INCLUDED_ALG_REMEZ_GENERAL_H
+#define INCLUDED_ALG_REMEZ_GENERAL_H
+
+#include <stddef.h>
+#include <Grid/GridStd.h>
+
+#ifdef HAVE_LIBGMP
+#include "bigfloat.h"
+#else
+#include "bigfloat_double.h"
+#endif
+
+
+class AlgRemezGeneral{
+ public:
+  enum PolyType { Even, Odd, Full };
+
+ private:
+
+  // In GSL-style, pass the function as a function pointer. Any data required to evaluate the function is passed in as a void pointer
+  bigfloat (*f)(bigfloat x, void *data);
+  void *data;
+
+  // The approximation parameters
+  std::vector<bigfloat> param;
+  bigfloat norm;
+
+  // The number of non-zero terms in the numerator and denominator
+  int n, d;
+  // The numerator and denominator degree (i.e.  the largest power)
+  int pow_n, pow_d;
+  
+  // Specify if the numerator and/or denominator are odd/even polynomials
+  PolyType num_type;
+  PolyType den_type;
+  std::vector<int> num_pows; //contains the mapping, with -1 if not present
+  std::vector<int> den_pows;
+
+  // The bounds of the approximation
+  bigfloat apstrt, apwidt, apend;
+
+  // Variables used to calculate the approximation
+  int nd1, iter;
+  std::vector<bigfloat> xx;
+  std::vector<bigfloat> mm;
+  std::vector<bigfloat> step;
+
+  bigfloat delta, spread;
+  
+  // Variables used in search
+  std::vector<bigfloat> yy;
+
+  // Variables used in solving linear equations
+  std::vector<bigfloat> A;
+  std::vector<bigfloat> B;
+  std::vector<int> IPS;
+
+  // The number of equations we must solve at each iteration (n+d+1)
+  int neq;
+
+  // The precision of the GNU MP library
+  long prec;
+
+  // Initialize member variables associated with the polynomial's properties
+  void setupPolyProperties(int num_degree, int den_degree, PolyType num_type_in, PolyType den_type_in);
+
+  // Initial values of maximal and minmal errors
+  void initialGuess();
+
+  // Initialise step sizes
+  void stpini();
+
+  // Initialize the algorithm
+  void reinitializeAlgorithm();
+
+  // Solve the equations
+  void equations();
+
+  // Search for error maxima and minima
+  void search(); 
+
+  // Calculate function required for the approximation
+  inline bigfloat func(bigfloat x) const{
+    return f(x, data);
+  }
+
+  // Compute size and sign of the approximation error at x
+  bigfloat getErr(bigfloat x, int *sign) const;
+
+  // Solve the system AX=B   where X = param
+  int simq();
+
+  // Evaluate the rational form P(x)/Q(x) using coefficients from the solution vector param
+  bigfloat approx(bigfloat x) const;
+
+ public:
+  
+  AlgRemezGeneral(double lower, double upper, long prec,
+		  bigfloat (*f)(bigfloat x, void *data), void *data);
+
+  inline int getDegree(void) const{ 
+    assert(n==d);
+    return n;
+  }
+  // Reset the bounds of the approximation
+  inline void setBounds(double lower, double upper) {
+    apstrt = lower;
+    apend = upper;
+    apwidt = apend - apstrt;
+  }
+
+  // Get the bounds of the approximation
+  inline void getBounds(double &lower, double &upper) const{ 
+    lower=(double)apstrt;
+    upper=(double)apend;
+  }
+
+  // Run the algorithm to generate the rational approximation
+  double generateApprox(int num_degree, int den_degree, 
+			PolyType num_type, PolyType den_type,
+			const double tolerance = 1e-15, const int report_freq = 1000);
+  
+  inline double generateApprox(int num_degree, int den_degree, 
+			       const double tolerance = 1e-15, const int report_freq = 1000){
+    return generateApprox(num_degree, den_degree, Full, Full, tolerance, report_freq);
+  }
+  
+  // Evaluate the rational form P(x)/Q(x) using coefficients from the
+  // solution vector param
+  inline double evaluateApprox(double x) const{
+    return (double)approx((bigfloat)x);
+  }
+
+  // Evaluate the rational form Q(x)/P(x) using coefficients from the solution vector param
+  inline double evaluateInverseApprox(double x) const{
+    return 1.0/(double)approx((bigfloat)x);
+  }  
+
+  // Calculate function required for the approximation
+  inline double evaluateFunc(double x) const{
+    return (double)func((bigfloat)x);
+  }
+
+  // Calculate inverse function required for the approximation
+  inline double evaluateInverseFunc(double x) const{
+    return 1.0/(double)func((bigfloat)x);
+  }
+
+  // Dump csv of function, approx and error
+  void csv(std::ostream &os = std::cout) const;
+
+  // Get the coefficient of the term x^i in the numerator
+  inline double getCoeffNum(const int i) const{    
+    return num_pows[i] == -1 ? 0. : double(param[num_pows[i]]);
+  }
+  // Get the coefficient of the term x^i in the denominator
+  inline double getCoeffDen(const int i) const{ 
+    if(i == pow_d) return 1.0;
+    else return den_pows[i] == -1 ? 0. : double(param[den_pows[i]+n+1]); 
+  }
+};
+
+#endif

Grid/algorithms/approx/ZMobius.cc

@@ -0,0 +1,183 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/approx/ZMobius.cc
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/algorithms/approx/ZMobius.h>
+#include <Grid/algorithms/approx/RemezGeneral.h>
+
+NAMESPACE_BEGIN(Grid);
+NAMESPACE_BEGIN(Approx);
+
+//Compute the tanh approximation
+inline double epsilonMobius(const double x, const std::vector<ComplexD> &w){
+  int Ls = w.size();
+
+  ComplexD fxp = 1., fmp = 1.;
+  for(int i=0;i<Ls;i++){
+    fxp = fxp * ( w[i] + x );
+    fmp = fmp * ( w[i] - x );
+  }
+  return ((fxp - fmp)/(fxp + fmp)).real();
+}
+inline double epsilonMobius(const double x, const std::vector<RealD> &w){
+  int Ls = w.size();
+
+  double fxp = 1., fmp = 1.;
+  for(int i=0;i<Ls;i++){
+    fxp = fxp * ( w[i] + x );
+    fmp = fmp * ( w[i] - x );
+  }
+  return (fxp - fmp)/(fxp + fmp);
+}
+
+
+
+//Compute the tanh approximation in a form suitable for the Remez
+bigfloat epsilonMobius(bigfloat x, void* data){
+  const std::vector<RealD> &omega = *( (std::vector<RealD> const*)data );
+  bigfloat fxp(1.0);
+  bigfloat fmp(1.0);
+
+  for(int i=0;i<omega.size();i++){
+    fxp = fxp * ( bigfloat(omega[i]) + x);
+    fmp = fmp * ( bigfloat(omega[i]) - x);
+  }
+  return (fxp - fmp)/(fxp + fmp);
+}
+
+//Compute the Zmobius Omega parameters suitable for eigenvalue range   -lambda_bound <= lambda <= lambda_bound
+//Note omega_i = 1/(b_i + c_i)   where b_i and c_i are the Mobius parameters
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out,
+			 const std::vector<RealD> &omega_in, const int Ls_in,
+			 const RealD lambda_bound){
+  assert(omega_in.size() == Ls_in);
+  omega_out.resize(Ls_out);
+
+  //Use the Remez algorithm to generate the appropriate rational polynomial
+  //For odd polynomial, to satisfy Haar condition must take either positive or negative half of range (cf https://arxiv.org/pdf/0803.0439.pdf page 6)  
+  AlgRemezGeneral remez(0, lambda_bound, 64, &epsilonMobius, (void*)&omega_in); 
+  remez.generateApprox(Ls_out-1, Ls_out,AlgRemezGeneral::Odd, AlgRemezGeneral::Even, 1e-15, 100);
+  remez.csv(std::cout);
+
+  //The rational approximation has the form  [ f(x) - f(-x) ] / [ f(x) + f(-x) ]  where  f(x) = \Prod_{i=0}^{L_s-1} ( \omega_i + x )
+  //cf https://academiccommons.columbia.edu/doi/10.7916/D8T72HD7  pg 102
+  //omega_i are therefore the negative of the complex roots of f(x)
+
+  //We can find the roots by recognizing that the eigenvalues of a matrix A are the roots of the characteristic polynomial
+  // \rho(\lambda) = det( A - \lambda I )    where I is the unit matrix
+  //The matrix whose characteristic polynomial is an arbitrary monic polynomial a0 + a1 x + a2 x^2 + ... x^n   is the companion matrix 
+  // A = | 0    1   0    0 0 .... 0 |
+  //     | 0    0   1    0 0 .... 0 |
+  //     | :    :   :    : :      : |
+  //     | 0    0   0    0 0      1
+  //     | -a0 -a1 -a2  ...  ... -an|
+
+
+  //Note the Remez defines the largest power to have unit coefficient
+  std::vector<RealD> coeffs(Ls_out+1);
+  for(int i=0;i<Ls_out+1;i+=2) coeffs[i] = coeffs[i] = remez.getCoeffDen(i); //even powers
+  for(int i=1;i<Ls_out+1;i+=2) coeffs[i] = coeffs[i] = remez.getCoeffNum(i); //odd powers
+
+  std::vector<std::complex<RealD> > roots(Ls_out);
+
+  //Form the companion matrix
+  Eigen::MatrixXd compn(Ls_out,Ls_out);
+  for(int i=0;i<Ls_out-1;i++) compn(i,0) = 0.;
+  compn(Ls_out - 1, 0) = -coeffs[0];
+  
+  for(int j=1;j<Ls_out;j++){
+    for(int i=0;i<Ls_out-1;i++) compn(i,j) = i == j-1 ? 1. : 0.;
+    compn(Ls_out - 1, j) = -coeffs[j];
+  }
+
+  //Eigensolve
+  Eigen::EigenSolver<Eigen::MatrixXd> slv(compn, false);
+
+  const auto & ev = slv.eigenvalues();
+  for(int i=0;i<Ls_out;i++)
+    omega_out[i] = -ev(i);
+
+  //Sort ascending (smallest at start of vector!)
+  std::sort(omega_out.begin(), omega_out.end(), 
+	    [&](const ComplexD &a, const ComplexD &b){ return a.real() < b.real() || (a.real() == b.real() && a.imag() < b.imag()); });
+
+  //McGlynn thesis pg 122 suggest improved iteration counts if magnitude of omega diminishes towards the center of the 5th dimension
+  std::vector<ComplexD> omega_tmp = omega_out;
+  int s_low=0, s_high=Ls_out-1, ss=0;
+  for(int s_from = Ls_out-1; s_from >= 0; s_from--){ //loop from largest omega
+    int s_to;
+    if(ss % 2 == 0){
+      s_to = s_low++;
+    }else{
+      s_to = s_high--;
+    }
+    omega_out[s_to] = omega_tmp[s_from];
+    ++ss;
+  }
+  
+  std::cout << "Resulting omega_i:" << std::endl;  
+  for(int i=0;i<Ls_out;i++)
+    std::cout << omega_out[i] << std::endl;
+
+  std::cout << "Test result matches the approximate polynomial found by the Remez" << std::endl;
+  std::cout << "<x> <remez approx> <poly approx> <diff poly approx remez approx> <exact> <diff poly approx exact>\n";
+  
+  int npt = 60;
+  double dlt = lambda_bound/double(npt-1);
+
+  for (int i =0; i<npt; i++){
+    double x = i*dlt;
+    double r = remez.evaluateApprox(x);
+    double p = epsilonMobius(x, omega_out);
+    double e = epsilonMobius(x, omega_in);
+
+    std::cout << x<< " " << r << " " << p <<" " <<r-p << " " << e << " " << e-p << std::endl;
+  }
+
+}
+  
+//mobius_param = b+c   with b-c=1
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound){
+  std::vector<RealD> omega_in(Ls_in, 1./mobius_param);
+  computeZmobiusOmega(omega_out, Ls_out, omega_in, Ls_in, lambda_bound);
+}
+
+//ZMobius class takes  gamma_i = (b+c) omega_i as its input, where b, c are factored out
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, 
+			 const RealD mobius_param_out, const int Ls_out, 
+			 const RealD mobius_param_in, const int Ls_in,
+			 const RealD lambda_bound){
+  computeZmobiusOmega(gamma_out, Ls_out, mobius_param_in, Ls_in, lambda_bound);
+  for(int i=0;i<Ls_out;i++) gamma_out[i] = gamma_out[i] * mobius_param_out;
+}
+//Assumes mobius_param_out == mobius_param_in
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound){
+  computeZmobiusGamma(gamma_out, mobius_param, Ls_out, mobius_param, Ls_in, lambda_bound);
+}
+
+NAMESPACE_END(Approx);
+NAMESPACE_END(Grid);

Grid/algorithms/approx/ZMobius.h

@@ -0,0 +1,57 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/approx/ZMobius.h
+
+    Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef GRID_ZMOBIUS_APPROX_H
+#define GRID_ZMOBIUS_APPROX_H
+
+#include <Grid/GridCore.h>
+
+NAMESPACE_BEGIN(Grid);
+NAMESPACE_BEGIN(Approx);
+
+//Compute the Zmobius Omega parameters suitable for eigenvalue range   -lambda_bound <= lambda <= lambda_bound
+//Note omega_i = 1/(b_i + c_i)   where b_i and c_i are the Mobius parameters
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out,
+			 const std::vector<RealD> &omega_in, const int Ls_in,
+			 const RealD lambda_bound);
+  
+//mobius_param = b+c   with b-c=1
+void computeZmobiusOmega(std::vector<ComplexD> &omega_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound);
+
+//ZMobius class takes  gamma_i = (b+c) omega_i as its input, where b, c are factored out
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, 
+			 const RealD mobius_param_out, const int Ls_out, 
+			 const RealD mobius_param_in, const int Ls_in,
+			 const RealD lambda_bound);
+
+//Assumes mobius_param_out == mobius_param_in
+void computeZmobiusGamma(std::vector<ComplexD> &gamma_out, const int Ls_out, const RealD mobius_param, const int Ls_in, const RealD lambda_bound);
+
+NAMESPACE_END(Approx);
+NAMESPACE_END(Grid);
+
+#endif

Grid/algorithms/approx/bigfloat_double.h

@@ -25,6 +25,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
+
+#ifndef INCLUDED_BIGFLOAT_DOUBLE_H
+#define INCLUDED_BIGFLOAT_DOUBLE_H
+
 #include <math.h>
 
 typedef double mfloat; 
@@ -186,4 +190,6 @@ public:
   //  friend bigfloat& random(void);
 };
 
+#endif
+
 

Grid/algorithms/iterative/BiCGSTAB.h

@@ -0,0 +1,222 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/algorithms/iterative/BiCGSTAB.h
+
+Copyright (C) 2015
+
+Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+Author: juettner <juettner@soton.ac.uk>
+Author: David Murphy <djmurphy@mit.edu>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef GRID_BICGSTAB_H
+#define GRID_BICGSTAB_H
+
+NAMESPACE_BEGIN(Grid);
+
+/////////////////////////////////////////////////////////////
+// Base classes for iterative processes based on operators
+// single input vec, single output vec.
+/////////////////////////////////////////////////////////////
+
+template <class Field>
+class BiCGSTAB : public OperatorFunction<Field> 
+{
+  public:
+    using OperatorFunction<Field>::operator();
+    
+    bool ErrorOnNoConverge;  // throw an assert when the CG fails to converge.
+                             // Defaults true.
+    RealD Tolerance;
+    Integer MaxIterations;
+    Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  
+    BiCGSTAB(RealD tol, Integer maxit, bool err_on_no_conv = true) : 
+      Tolerance(tol), MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv){};
+
+    void operator()(LinearOperatorBase<Field>& Linop, const Field& src, Field& psi) 
+    {
+      psi.Checkerboard() = src.Checkerboard();
+      conformable(psi, src);
+
+      RealD cp(0), rho(1), rho_prev(0), alpha(1), beta(0), omega(1);
+      RealD a(0), bo(0), b(0), ssq(0);
+
+      Field p(src);
+      Field r(src);
+      Field rhat(src);
+      Field v(src);
+      Field s(src);
+      Field t(src);
+      Field h(src);
+
+      v = Zero();
+      p = Zero();
+
+      // Initial residual computation & set up
+      RealD guess = norm2(psi);
+      assert(std::isnan(guess) == 0);
+    
+      Linop.Op(psi, v);
+      b = norm2(v);
+
+      r = src - v;
+      rhat = r;
+      a = norm2(r);
+      ssq = norm2(src);
+
+      std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: guess " << guess << std::endl;
+      std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB:   src " << ssq << std::endl;
+      std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB:    mp " << b << std::endl;
+      std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB:     r " << a << std::endl;
+
+      RealD rsq = Tolerance * Tolerance * ssq;
+
+      // Check if guess is really REALLY good :)
+      if(a <= rsq){ return; }
+
+      std::cout << GridLogIterative << std::setprecision(8) << "BiCGSTAB: k=0 residual " << a << " target " << rsq << std::endl;
+
+      GridStopWatch LinalgTimer;
+      GridStopWatch InnerTimer;
+      GridStopWatch AxpyNormTimer;
+      GridStopWatch LinearCombTimer;
+      GridStopWatch MatrixTimer;
+      GridStopWatch SolverTimer;
+
+      SolverTimer.Start();
+      int k;
+      for (k = 1; k <= MaxIterations; k++) 
+      {
+        rho_prev = rho;
+
+        LinalgTimer.Start();
+        InnerTimer.Start();
+        ComplexD Crho  = innerProduct(rhat,r);
+        InnerTimer.Stop();
+        rho = Crho.real();
+
+        beta = (rho / rho_prev) * (alpha / omega);
+
+        LinearCombTimer.Start();
+        bo = beta * omega;
+        auto p_v = p.View();
+        auto r_v = r.View();
+        auto v_v = v.View();
+        accelerator_for(ss, p_v.size(), Field::vector_object::Nsimd(),{
+          coalescedWrite(p_v[ss], beta*p_v(ss) - bo*v_v(ss) + r_v(ss));
+        });
+        LinearCombTimer.Stop();
+        LinalgTimer.Stop();
+
+        MatrixTimer.Start();
+        Linop.Op(p,v);
+        MatrixTimer.Stop();
+
+        LinalgTimer.Start();
+        InnerTimer.Start();
+        ComplexD Calpha = innerProduct(rhat,v);
+        InnerTimer.Stop();
+        alpha = rho / Calpha.real();
+
+        LinearCombTimer.Start();
+        auto h_v = h.View();
+        auto psi_v = psi.View();
+        accelerator_for(ss, h_v.size(), Field::vector_object::Nsimd(),{
+          coalescedWrite(h_v[ss], alpha*p_v(ss) + psi_v(ss));
+        });
+        
+        auto s_v = s.View();
+        accelerator_for(ss, s_v.size(), Field::vector_object::Nsimd(),{
+          coalescedWrite(s_v[ss], -alpha*v_v(ss) + r_v(ss));
+        });
+        LinearCombTimer.Stop();
+        LinalgTimer.Stop();
+
+        MatrixTimer.Start();
+        Linop.Op(s,t);
+        MatrixTimer.Stop();
+
+        LinalgTimer.Start();
+        InnerTimer.Start();
+        ComplexD Comega = innerProduct(t,s);
+        InnerTimer.Stop();
+        omega = Comega.real() / norm2(t);
+
+        LinearCombTimer.Start();
+        auto t_v = t.View();
+        accelerator_for(ss, psi_v.size(), Field::vector_object::Nsimd(),{
+          coalescedWrite(psi_v[ss], h_v(ss) + omega * s_v(ss));
+          coalescedWrite(r_v[ss], -omega * t_v(ss) + s_v(ss));
+        });
+        LinearCombTimer.Stop();
+
+        cp = norm2(r);
+        LinalgTimer.Stop();
+
+        std::cout << GridLogIterative << "BiCGSTAB: Iteration " << k << " residual " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
+
+        // Stopping condition
+        if(cp <= rsq) 
+        {
+          SolverTimer.Stop();
+          Linop.Op(psi, v);
+          p = v - src;
+
+          RealD srcnorm = sqrt(norm2(src));
+          RealD resnorm = sqrt(norm2(p));
+          RealD true_residual = resnorm / srcnorm;
+
+          std::cout << GridLogMessage << "BiCGSTAB Converged on iteration " << k << std::endl;
+          std::cout << GridLogMessage << "\tComputed residual " << sqrt(cp/ssq) << std::endl;
+          std::cout << GridLogMessage << "\tTrue residual " << true_residual << std::endl;
+          std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
+
+          std::cout << GridLogMessage << "Time breakdown " << std::endl;
+          std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() << std::endl;
+          std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() << std::endl;
+          std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() << std::endl;
+          std::cout << GridLogMessage << "\tInner      " << InnerTimer.Elapsed() << std::endl;
+          std::cout << GridLogMessage << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() << std::endl;
+          std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() << std::endl;
+
+          if(ErrorOnNoConverge){ assert(true_residual / Tolerance < 10000.0); }
+
+          IterationsToComplete = k;	
+
+          return;
+        }
+      }
+      
+      std::cout << GridLogMessage << "BiCGSTAB did NOT converge" << std::endl;
+
+      if(ErrorOnNoConverge){ assert(0); }
+      IterationsToComplete = k;
+    }
+};
+
+NAMESPACE_END(Grid);
+
+#endif

Grid/algorithms/iterative/BiCGSTABMixedPrec.h

@@ -0,0 +1,158 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: ./lib/algorithms/iterative/BiCGSTABMixedPrec.h
+
+Copyright (C) 2015
+
+Author: Christopher Kelly <ckelly@phys.columbia.edu>
+Author: David Murphy <djmurphy@mit.edu>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef GRID_BICGSTAB_MIXED_PREC_H
+#define GRID_BICGSTAB_MIXED_PREC_H
+
+NAMESPACE_BEGIN(Grid);
+
+// Mixed precision restarted defect correction BiCGSTAB
+template<class FieldD, class FieldF, typename std::enable_if< getPrecision<FieldD>::value == 2, int>::type = 0, typename std::enable_if< getPrecision<FieldF>::value == 1, int>::type = 0> 
+class MixedPrecisionBiCGSTAB : public LinearFunction<FieldD> 
+{
+  public:                                                
+    RealD   Tolerance;
+    RealD   InnerTolerance; // Initial tolerance for inner CG. Defaults to Tolerance but can be changed
+    Integer MaxInnerIterations;
+    Integer MaxOuterIterations;
+    GridBase* SinglePrecGrid; // Grid for single-precision fields
+    RealD OuterLoopNormMult; // Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
+    LinearOperatorBase<FieldF> &Linop_f;
+    LinearOperatorBase<FieldD> &Linop_d;
+
+    Integer TotalInnerIterations; //Number of inner CG iterations
+    Integer TotalOuterIterations; //Number of restarts
+    Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
+
+    //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess
+    LinearFunction<FieldF> *guesser;
+    
+    MixedPrecisionBiCGSTAB(RealD tol, Integer maxinnerit, Integer maxouterit, GridBase* _sp_grid, 
+        LinearOperatorBase<FieldF>& _Linop_f, LinearOperatorBase<FieldD>& _Linop_d) : 
+      Linop_f(_Linop_f), Linop_d(_Linop_d), Tolerance(tol), InnerTolerance(tol), MaxInnerIterations(maxinnerit), 
+      MaxOuterIterations(maxouterit), SinglePrecGrid(_sp_grid), OuterLoopNormMult(100.), guesser(NULL) {};
+
+    void useGuesser(LinearFunction<FieldF>& g){
+      guesser = &g;
+    }
+  
+    void operator() (const FieldD& src_d_in, FieldD& sol_d)
+    {
+      TotalInnerIterations = 0;
+    
+      GridStopWatch TotalTimer;
+      TotalTimer.Start();
+      
+      int cb = src_d_in.Checkerboard();
+      sol_d.Checkerboard() = cb;
+      
+      RealD src_norm = norm2(src_d_in);
+      RealD stop = src_norm * Tolerance*Tolerance;
+
+      GridBase* DoublePrecGrid = src_d_in.Grid();
+      FieldD tmp_d(DoublePrecGrid);
+      tmp_d.Checkerboard() = cb;
+      
+      FieldD tmp2_d(DoublePrecGrid);
+      tmp2_d.Checkerboard() = cb;
+      
+      FieldD src_d(DoublePrecGrid);
+      src_d = src_d_in; //source for next inner iteration, computed from residual during operation
+      
+      RealD inner_tol = InnerTolerance;
+      
+      FieldF src_f(SinglePrecGrid);
+      src_f.Checkerboard() = cb;
+      
+      FieldF sol_f(SinglePrecGrid);
+      sol_f.Checkerboard() = cb;
+      
+      BiCGSTAB<FieldF> CG_f(inner_tol, MaxInnerIterations);
+      CG_f.ErrorOnNoConverge = false;
+
+      GridStopWatch InnerCGtimer;
+
+      GridStopWatch PrecChangeTimer;
+      
+      Integer &outer_iter = TotalOuterIterations; //so it will be equal to the final iteration count
+        
+      for(outer_iter = 0; outer_iter < MaxOuterIterations; outer_iter++)
+      {
+        // Compute double precision rsd and also new RHS vector.
+        Linop_d.Op(sol_d, tmp_d);
+        RealD norm = axpy_norm(src_d, -1., tmp_d, src_d_in); //src_d is residual vector
+        
+        std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration " << outer_iter << " residual " << norm << " target " << stop << std::endl;
+
+        if(norm < OuterLoopNormMult * stop){
+          std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Outer iteration converged on iteration " << outer_iter << std::endl;
+          break;
+        }
+        while(norm * inner_tol * inner_tol < stop){ inner_tol *= 2; } // inner_tol = sqrt(stop/norm) ??
+
+        PrecChangeTimer.Start();
+        precisionChange(src_f, src_d);
+        PrecChangeTimer.Stop();
+        
+        sol_f = Zero();
+
+        //Optionally improve inner solver guess (eg using known eigenvectors)
+        if(guesser != NULL){ (*guesser)(src_f, sol_f); }
+
+        //Inner CG
+        CG_f.Tolerance = inner_tol;
+        InnerCGtimer.Start();
+        CG_f(Linop_f, src_f, sol_f);
+        InnerCGtimer.Stop();
+        TotalInnerIterations += CG_f.IterationsToComplete;
+        
+        //Convert sol back to double and add to double prec solution
+        PrecChangeTimer.Start();
+        precisionChange(tmp_d, sol_f);
+        PrecChangeTimer.Stop();
+        
+        axpy(sol_d, 1.0, tmp_d, sol_d);
+      }
+      
+      //Final trial CG
+      std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Starting final patch-up double-precision solve" << std::endl;
+      
+      BiCGSTAB<FieldD> CG_d(Tolerance, MaxInnerIterations);
+      CG_d(Linop_d, src_d_in, sol_d);
+      TotalFinalStepIterations = CG_d.IterationsToComplete;
+
+      TotalTimer.Stop();
+      std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Inner CG iterations " << TotalInnerIterations << " Restarts " << TotalOuterIterations << " Final CG iterations " << TotalFinalStepIterations << std::endl;
+      std::cout << GridLogMessage << "MixedPrecisionBiCGSTAB: Total time " << TotalTimer.Elapsed() << " Precision change " << PrecChangeTimer.Elapsed() << " Inner CG total " << InnerCGtimer.Elapsed() << std::endl;
+  }
+};
+
+NAMESPACE_END(Grid);
+
+#endif

Grid/algorithms/iterative/BlockConjugateGradient.h

@@ -52,6 +52,7 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
   Integer PrintInterval; //GridLogMessages or Iterative
+  RealD TrueResidual;
   
   BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true)
     : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100)
@@ -306,7 +307,8 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
 
       Linop.HermOp(X, AD);
       AD = AD-B;
-      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(norm2(AD)/norm2(B)) <<std::endl;
+      TrueResidual = std::sqrt(norm2(AD)/norm2(B));
+      std::cout << GridLogMessage <<"\tTrue residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
@@ -442,7 +444,8 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
 
       Linop.HermOp(Psi, AP);
       AP = AP-Src;
-      std::cout <<GridLogMessage << "\tTrue residual is " << std::sqrt(norm2(AP)/norm2(Src)) <<std::endl;
+      TrueResidual = std::sqrt(norm2(AP)/norm2(Src));
+      std::cout <<GridLogMessage << "\tTrue residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;
@@ -653,7 +656,7 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
       if ( rr > max_resid ) max_resid = rr;
     }
 
-    std::cout << GridLogIterative << "\t Block Iteration "<<k<<" ave resid "<< sqrt(rrsum/sssum) << " max "<< sqrt(max_resid) <<std::endl;
+    std::cout << GridLogIterative << "\t Block Iteration "<<k<<" ave resid "<< std::sqrt(rrsum/sssum) << " max "<< std::sqrt(max_resid) <<std::endl;
 
     if ( max_resid < Tolerance*Tolerance ) { 
 
@@ -668,7 +671,8 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
 
       for(int b=0;b<Nblock;b++) Linop.HermOp(X[b], AD[b]);
       for(int b=0;b<Nblock;b++) AD[b] = AD[b]-B[b];
-      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(normv(AD)/normv(B)) <<std::endl;
+      TrueResidual = std::sqrt(normv(AD)/normv(B));
+      std::cout << GridLogMessage << "\tTrue residual is " << TrueResidual <<std::endl;
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;
       std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed()     <<std::endl;

Grid/algorithms/iterative/ConjugateGradient.h

@@ -49,6 +49,7 @@ public:
   RealD Tolerance;
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  RealD TrueResidual;
   
   ConjugateGradient(RealD tol, Integer maxit, bool err_on_no_conv = true)
     : Tolerance(tol),
@@ -81,6 +82,14 @@ public:
     cp = a;
     ssq = norm2(src);
 
+    // Handle trivial case of zero src
+    if (ssq == 0.){
+      psi = Zero();
+      IterationsToComplete = 1;
+      TrueResidual = 0.;
+      return;
+    }
+
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient: guess " << guess << std::endl;
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient:   src " << ssq << std::endl;
     std::cout << GridLogIterative << std::setprecision(8) << "ConjugateGradient:    mp " << d << std::endl;
@@ -92,6 +101,7 @@ public:
 
     // Check if guess is really REALLY good :)
     if (cp <= rsq) {
+      TrueResidual = std::sqrt(a/ssq);
       std::cout << GridLogMessage << "ConjugateGradient guess is converged already " << std::endl;
       IterationsToComplete = 0;	
       return;
@@ -141,7 +151,7 @@ public:
       LinalgTimer.Stop();
 
       std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k
-                << " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
+                << " residual " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
 
       // Stopping condition
       if (cp <= rsq) {
@@ -169,10 +179,17 @@ public:
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 
 	IterationsToComplete = k;	
+	TrueResidual = true_residual;
 
         return;
       }
     }
+    // Failed. Calculate true residual before giving up                                                         
+    Linop.HermOpAndNorm(psi, mmp, d, qq);
+    p = mmp - src;
+
+    TrueResidual = sqrt(norm2(p)/ssq);
+
     std::cout << GridLogMessage << "ConjugateGradient did NOT converge "<<k<<" / "<< MaxIterations<< std::endl;
 
     if (ErrorOnNoConverge) assert(0);

Grid/algorithms/iterative/ConjugateGradientMultiShift.h

@@ -46,15 +46,19 @@ public:
 
   RealD   Tolerance;
   Integer MaxIterations;
-    Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  std::vector<int> IterationsToCompleteShift;  // Iterations for this shift
   int verbose;
   MultiShiftFunction shifts;
+  std::vector<RealD> TrueResidualShift;
 
   ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : 
     MaxIterations(maxit),
     shifts(_shifts)
   { 
     verbose=1;
+    IterationsToCompleteShift.resize(_shifts.order);
+    TrueResidualShift.resize(_shifts.order);
   }
 
   void operator() (LinearOperatorBase<Field> &Linop, const Field &src, Field &psi)
@@ -125,6 +129,17 @@ public:
     // Residuals "r" are src
     // First search direction "p" is also src
     cp = norm2(src);
+
+    // Handle trivial case of zero src.
+    if( cp == 0. ){
+      for(int s=0;s<nshift;s++){
+	psi[s] = Zero();
+	IterationsToCompleteShift[s] = 1;
+	TrueResidualShift[s] = 0.;
+      }
+      return;
+    }
+
     for(int s=0;s<nshift;s++){
       rsq[s] = cp * mresidual[s] * mresidual[s];
       std::cout<<GridLogMessage<<"ConjugateGradientMultiShift: shift "<<s
@@ -270,6 +285,7 @@ public:
       for(int s=0;s<nshift;s++){
       
 	if ( (!converged[s]) ){
+	  IterationsToCompleteShift[s] = k;
 	
 	  RealD css  = c * z[s][iz]* z[s][iz];
 	
@@ -299,7 +315,8 @@ public:
 	  axpy(r,-alpha[s],src,tmp);
 	  RealD rn = norm2(r);
 	  RealD cn = norm2(src);
-	  std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<<std::sqrt(rn/cn)<<std::endl;
+	  TrueResidualShift[s] = std::sqrt(rn/cn);
+	  std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<< TrueResidualShift[s] <<std::endl;
 	}
 
       std::cout << GridLogMessage << "Time Breakdown "<<std::endl;

Grid/algorithms/iterative/SchurRedBlack.h

@@ -405,6 +405,70 @@ namespace Grid {
     }
   };
 
+  template<class Field> class NonHermitianSchurRedBlackDiagMooeeSolve : public SchurRedBlackBase<Field> 
+  {
+    public:
+      typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
+
+      NonHermitianSchurRedBlackDiagMooeeSolve(OperatorFunction<Field>& RBSolver, const bool initSubGuess = false,
+          const bool _solnAsInitGuess = false)  
+      : SchurRedBlackBase<Field>(RBSolver, initSubGuess, _solnAsInitGuess) {};
+
+      //////////////////////////////////////////////////////
+      // Override RedBlack specialisation
+      //////////////////////////////////////////////////////
+      virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o)
+      {
+        GridBase* grid  = _Matrix.RedBlackGrid();
+        GridBase* fgrid = _Matrix.Grid();
+
+        Field  tmp(grid);
+        Field Mtmp(grid);
+
+        pickCheckerboard(Even, src_e, src);
+        pickCheckerboard(Odd , src_o, src);
+
+        /////////////////////////////////////////////////////
+        // src_o = Mdag * (source_o - Moe MeeInv source_e)
+        /////////////////////////////////////////////////////
+        _Matrix.MooeeInv(src_e, tmp);   assert(   tmp.Checkerboard() == Even );
+        _Matrix.Meooe   (tmp, Mtmp);    assert(  Mtmp.Checkerboard() == Odd  );     
+        src_o -= Mtmp;                  assert( src_o.Checkerboard() == Odd  );     
+      }
+      
+      virtual void RedBlackSolution(Matrix& _Matrix, const Field& sol_o, const Field& src_e, Field& sol)
+      {
+        GridBase* grid  = _Matrix.RedBlackGrid();
+        GridBase* fgrid = _Matrix.Grid();
+
+        Field     tmp(grid);
+        Field   sol_e(grid);
+        Field src_e_i(grid);
+        
+        ///////////////////////////////////////////////////
+        // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+        ///////////////////////////////////////////////////
+        _Matrix.Meooe(sol_o, tmp);         assert(     tmp.Checkerboard() == Even );
+        src_e_i = src_e - tmp;             assert( src_e_i.Checkerboard() == Even );
+        _Matrix.MooeeInv(src_e_i, sol_e);  assert(   sol_e.Checkerboard() == Even );
+       
+        setCheckerboard(sol, sol_e); assert( sol_e.Checkerboard() == Even );
+        setCheckerboard(sol, sol_o); assert( sol_o.Checkerboard() == Odd  );
+      }
+
+      virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o)
+      {
+        NonHermitianSchurDiagMooeeOperator<Matrix,Field> _OpEO(_Matrix);
+        this->_HermitianRBSolver(_OpEO, src_o, sol_o);  assert(sol_o.Checkerboard() == Odd);
+      }
+
+      virtual void RedBlackSolve(Matrix& _Matrix, const std::vector<Field>& src_o, std::vector<Field>& sol_o)
+      {
+        NonHermitianSchurDiagMooeeOperator<Matrix,Field> _OpEO(_Matrix);
+        this->_HermitianRBSolver(_OpEO, src_o, sol_o); 
+      }
+  };
+
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
   // Site diagonal is identity, right preconditioned by Mee^inv
   // ( 1 - Meo Moo^inv Moe Mee^inv  ) phi =( 1 - Meo Moo^inv Moe Mee^inv  ) Mee psi =  = eta  = eta
@@ -482,5 +546,76 @@ namespace Grid {
       this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); 
     }
   };
+
+  template<class Field> class NonHermitianSchurRedBlackDiagTwoSolve : public SchurRedBlackBase<Field> 
+  {
+    public:
+      typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
+
+      /////////////////////////////////////////////////////
+      // Wrap the usual normal equations Schur trick
+      /////////////////////////////////////////////////////
+      NonHermitianSchurRedBlackDiagTwoSolve(OperatorFunction<Field>& RBSolver, const bool initSubGuess = false,
+          const bool _solnAsInitGuess = false)  
+      : SchurRedBlackBase<Field>(RBSolver, initSubGuess, _solnAsInitGuess) {};
+
+      virtual void RedBlackSource(Matrix& _Matrix, const Field& src, Field& src_e, Field& src_o)
+      {
+        GridBase* grid  = _Matrix.RedBlackGrid();
+        GridBase* fgrid = _Matrix.Grid();
+
+        Field  tmp(grid);
+        Field Mtmp(grid);
+
+        pickCheckerboard(Even, src_e, src);
+        pickCheckerboard(Odd , src_o, src);
+      
+        /////////////////////////////////////////////////////
+        // src_o = Mdag * (source_o - Moe MeeInv source_e)
+        /////////////////////////////////////////////////////
+        _Matrix.MooeeInv(src_e, tmp);   assert(   tmp.Checkerboard() == Even );
+        _Matrix.Meooe   (tmp, Mtmp);    assert(  Mtmp.Checkerboard() == Odd  );     
+        src_o -= Mtmp;                  assert( src_o.Checkerboard() == Odd  );     
+      }
+
+      virtual void RedBlackSolution(Matrix& _Matrix, const Field& sol_o, const Field& src_e, Field& sol)
+      {
+        GridBase* grid  = _Matrix.RedBlackGrid();
+        GridBase* fgrid = _Matrix.Grid();
+
+        Field sol_o_i(grid);
+        Field     tmp(grid);
+        Field   sol_e(grid);
+
+        ////////////////////////////////////////////////
+        // MooeeInv due to pecond
+        ////////////////////////////////////////////////
+        _Matrix.MooeeInv(sol_o, tmp);
+        sol_o_i = tmp;
+
+        ///////////////////////////////////////////////////
+        // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+        ///////////////////////////////////////////////////
+        _Matrix.Meooe(sol_o_i, tmp);    assert(   tmp.Checkerboard() == Even );
+        tmp = src_e - tmp;              assert( src_e.Checkerboard() == Even );
+        _Matrix.MooeeInv(tmp, sol_e);   assert( sol_e.Checkerboard() == Even );
+       
+        setCheckerboard(sol, sol_e);    assert(   sol_e.Checkerboard() == Even );
+        setCheckerboard(sol, sol_o_i);  assert( sol_o_i.Checkerboard() == Odd  );
+      };
+
+      virtual void RedBlackSolve(Matrix& _Matrix, const Field& src_o, Field& sol_o)
+      {
+        NonHermitianSchurDiagTwoOperator<Matrix,Field> _OpEO(_Matrix);
+        this->_HermitianRBSolver(_OpEO, src_o, sol_o);
+      };
+
+      virtual void RedBlackSolve(Matrix& _Matrix, const std::vector<Field>& src_o,  std::vector<Field>& sol_o)
+      {
+        NonHermitianSchurDiagTwoOperator<Matrix,Field> _OpEO(_Matrix);
+        this->_HermitianRBSolver(_OpEO, src_o, sol_o); 
+      }
+  };
 }
+
 #endif

Grid/lattice/Lattice_peekpoke.h

@@ -156,7 +156,7 @@ void peekSite(sobj &s,const Lattice<vobj> &l,const Coordinate &site){
 // Peek a scalar object from the SIMD array
 //////////////////////////////////////////////////////////
 template<class vobj,class sobj>
-accelerator_inline void peekLocalSite(sobj &s,const Lattice<vobj> &l,Coordinate &site){
+inline void peekLocalSite(sobj &s,const Lattice<vobj> &l,Coordinate &site){
         
   GridBase *grid = l.Grid();
 
@@ -185,7 +185,7 @@ accelerator_inline void peekLocalSite(sobj &s,const Lattice<vobj> &l,Coordinate
 };
 
 template<class vobj,class sobj>
-accelerator_inline void pokeLocalSite(const sobj &s,Lattice<vobj> &l,Coordinate &site){
+inline void pokeLocalSite(const sobj &s,Lattice<vobj> &l,Coordinate &site){
 
   GridBase *grid=l.Grid();
 

Grid/perfmon/PerfCount.h

@@ -95,7 +95,8 @@ inline uint64_t cyclecount(void){
 }
 #elif defined __x86_64__
 inline uint64_t cyclecount(void){ 
-  return __rdtsc();
+  uint64_t ret = __rdtsc();
+  return (uint64_t)ret;
 }
 #else
 

Grid/qcd/QCD.h

@@ -133,23 +133,23 @@ typedef iSpinColourMatrix<vComplex >    vSpinColourMatrix;
 typedef iSpinColourMatrix<vComplexF>    vSpinColourMatrixF;
 typedef iSpinColourMatrix<vComplexD>    vSpinColourMatrixD;
 
-    // SpinColourSpinColour matrix
+// SpinColourSpinColour matrix
-    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
-    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
-    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
 
-    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
-    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
-    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
+typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
 
-    // SpinColourSpinColour matrix
+// SpinColourSpinColour matrix
-    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
-    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
-    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
 
-    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
-    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
-    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
+typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
 
 // LorentzColour
 typedef iLorentzColourMatrix<Complex  > LorentzColourMatrix;
@@ -443,16 +443,16 @@ template<class vobj> void pokeLorentz(vobj &lhs,const decltype(peekIndex<Lorentz
 //////////////////////////////////////////////
 // Fermion <-> propagator assignements
 //////////////////////////////////////////////
-    //template <class Prop, class Ferm>
+//template <class Prop, class Ferm>
-    template <class Fimpl>
+template <class Fimpl>
-      void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c)
+void FermToProp(typename Fimpl::PropagatorField &p, const typename Fimpl::FermionField &f, const int s, const int c)
 {
   for(int j = 0; j < Ns; ++j)
     {
       auto pjs = peekSpin(p, j, s);
       auto fj  = peekSpin(f, j);
             
-            for(int i = 0; i < Fimpl::Dimension; ++i)
+      for(int i = 0; i < Fimpl::Dimension; ++i)
 	{
 	  pokeColour(pjs, peekColour(fj, i), i, c);
 	}
@@ -460,16 +460,16 @@ template<class vobj> void pokeLorentz(vobj &lhs,const decltype(peekIndex<Lorentz
     }
 }
     
-    //template <class Prop, class Ferm>
+//template <class Prop, class Ferm>
-    template <class Fimpl>
+template <class Fimpl>
-      void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c)
+void PropToFerm(typename Fimpl::FermionField &f, const typename Fimpl::PropagatorField &p, const int s, const int c)
 {
   for(int j = 0; j < Ns; ++j)
     {
       auto pjs = peekSpin(p, j, s);
       auto fj  = peekSpin(f, j);
             
-            for(int i = 0; i < Fimpl::Dimension; ++i)
+      for(int i = 0; i < Fimpl::Dimension; ++i)
 	{
 	  pokeColour(fj, peekColour(pjs, i, c), i);
 	}

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -141,7 +141,33 @@ public:
   Vector<iSinglet<Simd> >  MatpInvDag;
   Vector<iSinglet<Simd> >  MatmInvDag;
 
+  ///////////////////////////////////////////////////////////////
+  // Conserved current utilities
+  ///////////////////////////////////////////////////////////////
+
+  // Virtual can't template
+  void ContractConservedCurrent(PropagatorField &q_in_1,
+				PropagatorField &q_in_2,
+				PropagatorField &q_out,
+				PropagatorField &phys_src,
+				Current curr_type, 
+				unsigned int mu);
+
+   void SeqConservedCurrent(PropagatorField &q_in,
+			   PropagatorField &q_out,
+			   PropagatorField &phys_src,
+			   Current curr_type,
+			   unsigned int mu,
+			   unsigned int tmin,
+			   unsigned int tmax,
+			   ComplexField &lattice_cmplx);
+
+  void ContractJ5q(PropagatorField &q_in,ComplexField &J5q);
+  void ContractJ5q(FermionField &q_in,ComplexField &J5q);
+
+  ///////////////////////////////////////////////////////////////
   // Constructors
+  ///////////////////////////////////////////////////////////////
   CayleyFermion5D(GaugeField &_Umu,
 		  GridCartesian         &FiveDimGrid,
 		  GridRedBlackCartesian &FiveDimRedBlackGrid,

Grid/qcd/action/fermion/FermionOperator.h

@@ -148,15 +148,19 @@ public:
   virtual void ContractConservedCurrent(PropagatorField &q_in_1,
 					PropagatorField &q_in_2,
 					PropagatorField &q_out,
+					PropagatorField &phys_src,
 					Current curr_type,
-					unsigned int mu)=0;
+					unsigned int mu)
+  {assert(0);};
   virtual void SeqConservedCurrent(PropagatorField &q_in, 
 				   PropagatorField &q_out,
+				   PropagatorField &phys_src,
 				   Current curr_type,
 				   unsigned int mu,
 				   unsigned int tmin, 
 				   unsigned int tmax,
-				   ComplexField &lattice_cmplx)=0;
+				   ComplexField &lattice_cmplx)
+  {assert(0);};
 
       // Only reimplemented in Wilson5D 
       // Default to just a zero correlation function

Grid/qcd/action/fermion/GparityWilsonImpl.h

@@ -38,6 +38,7 @@ public:
  static const bool isFundamental = Representation::isFundamental;
  static const int Nhcs = Options::Nhcs;
  static const bool LsVectorised=false;
+ static const bool isGparity=true;
 
  typedef ConjugateGaugeImpl< GaugeImplTypes<S,Dimension> > Gimpl;
  INHERIT_GIMPL_TYPES(Gimpl);
@@ -46,7 +47,7 @@ public:
  typedef typename Options::template PrecisionMapper<Simd>::LowerPrecVector SimdL;
       
  template <typename vtype> using iImplSpinor            = iVector<iVector<iVector<vtype, Dimension>, Ns>,   Ngp>;
- template <typename vtype> using iImplPropagator        = iVector<iMatrix<iMatrix<vtype, Dimension>, Ns>,   Ngp>;
+ template <typename vtype> using iImplPropagator        = iMatrix<iMatrix<iMatrix<vtype, Dimension>, Ns>,   Ngp>;
  template <typename vtype> using iImplHalfSpinor        = iVector<iVector<iVector<vtype, Dimension>, Nhs>,  Ngp>;
  template <typename vtype> using iImplHalfCommSpinor    = iVector<iVector<iVector<vtype, Dimension>, Nhcs>, Ngp>;
  template <typename vtype> using iImplDoubledGaugeField = iVector<iVector<iScalar<iMatrix<vtype, Dimension> >, Nds>, Ngp>;
@@ -80,6 +81,7 @@ public:
   {
     assert(0);
   } 
+
   template<class _Spinor>
   static accelerator_inline void multLink(_Spinor &phi, 
 					  const SiteDoubledGaugeField &U,
@@ -191,6 +193,16 @@ public:
 #endif   
   }
 
+
+  template<class _SpinorField>
+  inline void multLinkField(_SpinorField & out,
+			    const DoubledGaugeField &Umu,
+			    const _SpinorField & phi,
+			    int mu)
+  {
+    assert(0);
+  }
+
   template <class ref>
   static accelerator_inline void loadLinkElement(Simd &reg, ref &memory) 
   {

Grid/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -185,10 +185,12 @@ public:
   void ContractConservedCurrent(PropagatorField &q_in_1,
                                 PropagatorField &q_in_2,
                                 PropagatorField &q_out,
+                                PropagatorField &src,
                                 Current curr_type,
                                 unsigned int mu);
   void SeqConservedCurrent(PropagatorField &q_in,
                            PropagatorField &q_out,
+                           PropagatorField &srct,
                            Current curr_type,
                            unsigned int mu, 
                            unsigned int tmin,

Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -217,15 +217,17 @@ public:
   void ContractConservedCurrent(PropagatorField &q_in_1,
 				PropagatorField &q_in_2,
 				PropagatorField &q_out,
+				PropagatorField &src,
 				Current curr_type,
 				unsigned int mu);
   void SeqConservedCurrent(PropagatorField &q_in,
 			   PropagatorField &q_out,
+			   PropagatorField &src,
 			   Current curr_type,
 			   unsigned int mu, 
 			   unsigned int tmin,
-                             unsigned int tmax,
+			   unsigned int tmax,
-                 	     ComplexField &lattice_cmplx);
+			   ComplexField &lattice_cmplx);
 };
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/MADWF.h

@@ -40,6 +40,11 @@ inline void convert(const Fieldi &from,Fieldo &to)
   to=from;
 }
 
+struct MADWFinnerIterCallbackBase{
+  virtual void operator()(const RealD current_resid){}
+  virtual ~MADWFinnerIterCallbackBase(){}
+};
+
 template<class Matrixo,class Matrixi,class PVinverter,class SchurSolver, class Guesser> 
 class MADWF 
 {
@@ -56,24 +61,30 @@ class MADWF
 
   RealD target_resid;
   int   maxiter;
- public:
 
+  //operator() is called on "callback" at the end of every inner iteration. This allows for example the adjustment of the inner
+  //tolerance to speed up subsequent iteration
+  MADWFinnerIterCallbackBase* callback;
+  
+ public:
   MADWF(Matrixo &_Mato,
-	Matrixi &_Mati, 
+	Matrixi &_Mati,
-	PVinverter &_PauliVillarsSolvero, 
+	PVinverter &_PauliVillarsSolvero,
 	SchurSolver &_SchurSolveri,
 	Guesser & _Guesseri,
 	RealD resid,
-	int _maxiter) :
+	int _maxiter,
+	MADWFinnerIterCallbackBase* _callback = NULL) :
 
   Mato(_Mato),Mati(_Mati),
     SchurSolveri(_SchurSolveri),
-    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri)
+    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri),
-  {   
+    callback(_callback)
-    target_resid=resid;
+    {
-    maxiter     =_maxiter; 
+      target_resid=resid;
-  };
+      maxiter     =_maxiter;
-
+    };
+   
   void operator() (const FermionFieldo &src4,FermionFieldo &sol5)
   {
     std::cout << GridLogMessage<< " ************************************************" << std::endl;
@@ -177,6 +188,8 @@ class MADWF
        std::cout << GridLogMessage << "Residual " << i << ": " << resid  << std::endl;
        std::cout << GridLogMessage << "***************************************" <<std::endl;
 
+       if(callback != NULL) (*callback)(resid);       
+       
        if (resid < target_resid) {
 	 return;
        }

Grid/qcd/action/fermion/WilsonFermion.h

@@ -179,15 +179,17 @@ public:
   void ContractConservedCurrent(PropagatorField &q_in_1,
                                 PropagatorField &q_in_2,
                                 PropagatorField &q_out,
+                                PropagatorField &phys_src,
                                 Current curr_type,
                                 unsigned int mu);
   void SeqConservedCurrent(PropagatorField &q_in,
                            PropagatorField &q_out,
+                           PropagatorField &phys_src,
                            Current curr_type,
                            unsigned int mu, 
                            unsigned int tmin,
-                             unsigned int tmax,
+			   unsigned int tmax,
-			     ComplexField &lattice_cmplx);
+			   ComplexField &lattice_cmplx);
 };
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;

Grid/qcd/action/fermion/WilsonFermion5D.h

@@ -217,25 +217,7 @@ public:
     
   // Comms buffer
   std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  comm_buf;
-    
-  ///////////////////////////////////////////////////////////////
-  // Conserved current utilities
-  ///////////////////////////////////////////////////////////////
-  void ContractConservedCurrent(PropagatorField &q_in_1,
-				PropagatorField &q_in_2,
-				PropagatorField &q_out,
-				Current curr_type, 
-				unsigned int mu);
-  void SeqConservedCurrent(PropagatorField &q_in,
-			   PropagatorField &q_out,
-			   Current curr_type,
-			   unsigned int mu,
-			   unsigned int tmin,
-			   unsigned int tmax,
-			   ComplexField &lattice_cmplx);
 
-  void ContractJ5q(PropagatorField &q_in,ComplexField &J5q);
-  void ContractJ5q(FermionField &q_in,ComplexField &J5q);
 
 };
 

Grid/qcd/action/fermion/WilsonImpl.h

@@ -41,6 +41,7 @@ public:
   static const int Dimension = Representation::Dimension;
   static const bool isFundamental = Representation::isFundamental;
   static const bool LsVectorised=false;
+  static const bool isGparity=false;
   static const int Nhcs = Options::Nhcs;
 
   typedef PeriodicGaugeImpl<GaugeImplTypes<S, Dimension > > Gimpl;
@@ -98,8 +99,21 @@ public:
   {
     multLink(phi,U,chi,mu);
   }
-    
+
-      
+  template<class _SpinorField> 
+  inline void multLinkField(_SpinorField & out,
+			    const DoubledGaugeField &Umu,
+			    const _SpinorField & phi,
+			    int mu)
+  {
+    auto out_v= out.View();
+    auto phi_v= phi.View();
+    auto Umu_v= Umu.View();
+    thread_for(sss,out.Grid()->oSites(),{
+	multLink(out_v[sss],Umu_v[sss],phi_v[sss],mu);
+    });
+  }
+					   
   template <class ref>
   static accelerator_inline void loadLinkElement(Simd &reg, ref &memory) 
   {

Grid/qcd/action/fermion/WilsonKernels.h

@@ -66,41 +66,6 @@ public:
   static void DhopDirKernel(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf,
 			    int Ls, int Nsite, const FermionField &in, FermionField &out, int dirdisp, int gamma);
 
-  //////////////////////////////////////////////////////////////////////////////
-  // Utilities for inserting Wilson conserved current.
-  //////////////////////////////////////////////////////////////////////////////
-  static void ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1,
-                                       const SitePropagator &q_in_2,
-                                       SitePropagator &q_out,
-                                       DoubledGaugeFieldView &U,
-                                       unsigned int sU,
-                                       unsigned int mu,
-                                       bool switch_sign = false);
-
-  static void ContractConservedCurrentSiteBwd(const SitePropagator &q_in_1,
-                                       const SitePropagator &q_in_2,
-                                       SitePropagator &q_out,
-                                       DoubledGaugeFieldView &U,
-                                       unsigned int sU,
-                                       unsigned int mu,
-                                       bool switch_sign = false);
-
-  static void SeqConservedCurrentSiteFwd(const SitePropagator &q_in, 
-                                  SitePropagator &q_out,
-                                  DoubledGaugeFieldView &U,
-                                  unsigned int sU,
-                                  unsigned int mu,
-                                  vPredicate t_mask,
-                                  bool switch_sign = false);
-
-  static void SeqConservedCurrentSiteBwd(const SitePropagator &q_in,
-                                  SitePropagator &q_out,
-                                  DoubledGaugeFieldView &U,
-                                  unsigned int sU,
-                                  unsigned int mu,
-                                  vPredicate t_mask,
-                                  bool switch_sign = false);
-
 private:
 
   static accelerator_inline void DhopDirK(StencilView &st, DoubledGaugeFieldView &U,SiteHalfSpinor * buf,

Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h

@@ -588,6 +588,355 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,Vector<Coeff_t
   //  this->MooeeInternalCompute(1,inv,MatpInvDag,MatmInvDag);
 }
 
+
+template <class Impl>
+void CayleyFermion5D<Impl>::ContractJ5q(FermionField &q_in,ComplexField &J5q)
+{
+  conformable(this->GaugeGrid(), J5q.Grid());
+  conformable(q_in.Grid(), this->FermionGrid());
+  Gamma G5(Gamma::Algebra::Gamma5);
+  // 4d field
+  int Ls = this->Ls;
+  FermionField psi(this->GaugeGrid());
+  FermionField p_plus (this->GaugeGrid());
+  FermionField p_minus(this->GaugeGrid());
+  FermionField p(this->GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2-1 , 0);
+  ExtractSlice(p_minus, q_in, Ls/2   , 0);
+  p_plus = p_plus + G5*p_plus;
+  p_minus= p_minus - G5*p_minus;
+  p=0.5*(p_plus+p_minus);
+  J5q = localInnerProduct(p,p);
+}
+
+template <class Impl>
+void CayleyFermion5D<Impl>::ContractJ5q(PropagatorField &q_in,ComplexField &J5q)
+{
+  conformable(this->GaugeGrid(), J5q.Grid());
+  conformable(q_in.Grid(), this->FermionGrid());
+  Gamma G5(Gamma::Algebra::Gamma5);
+  // 4d field
+  int Ls = this->Ls;
+  PropagatorField psi(this->GaugeGrid());
+  PropagatorField p_plus (this->GaugeGrid());
+  PropagatorField p_minus(this->GaugeGrid());
+  PropagatorField p(this->GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2-1 , 0);
+  ExtractSlice(p_minus, q_in, Ls/2   , 0);
+  p_plus = p_plus + G5*p_plus;
+  p_minus= p_minus - G5*p_minus;
+  p=0.5*(p_plus+p_minus);
+  J5q = localInnerProduct(p,p);
+}
+
+#define Pp(Q) (0.5*(Q+g5*Q))
+#define Pm(Q) (0.5*(Q-g5*Q))
+#define Q_4d(Q) (Pm((Q)[0]) + Pp((Q)[Ls-1]))
+#define TopRowWithSource(Q) (phys_src + (1.0-mass)*Q_4d(Q))
+
+template <class Impl> 
+void CayleyFermion5D<Impl>::ContractConservedCurrent( PropagatorField &q_in_1,
+						      PropagatorField &q_in_2,
+						      PropagatorField &q_out,
+						      PropagatorField &phys_src,
+						      Current curr_type,
+						      unsigned int mu)
+{
+#ifndef GRID_NVCC
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT,
+    Gamma::Algebra::Gamma5
+  };
+
+  auto UGrid= this->GaugeGrid();
+  auto FGrid= this->FermionGrid();
+  RealD sgn=1.0;
+  if ( curr_type == Current::Axial ) sgn = -1.0;
+
+  int Ls = this->Ls;
+
+  std::vector<PropagatorField> L_Q(Ls,UGrid); 
+  std::vector<PropagatorField> R_Q(Ls,UGrid); 
+  for(int s=0;s<Ls;s++){
+    ExtractSlice(L_Q[s], q_in_1, s , 0);
+    ExtractSlice(R_Q[s], q_in_2, s , 0);
+  }
+
+  Gamma g5(Gamma::Algebra::Gamma5);
+  PropagatorField C(UGrid); 
+  PropagatorField p5d(UGrid); 
+  PropagatorField us_p5d(UGrid); 
+  PropagatorField gp5d(UGrid); 
+  PropagatorField gus_p5d(UGrid); 
+
+  PropagatorField L_TmLsGq0(UGrid); 
+  PropagatorField L_TmLsTmp(UGrid);
+  PropagatorField R_TmLsGq0(UGrid); 
+  PropagatorField R_TmLsTmp(UGrid);
+  {
+    PropagatorField TermA(UGrid);
+    PropagatorField TermB(UGrid);
+    PropagatorField TermC(UGrid);
+    PropagatorField TermD(UGrid);
+    TermA = (Pp(Q_4d(L_Q)));
+    TermB = (Pm(Q_4d(L_Q)));
+    TermC = (Pm(TopRowWithSource(L_Q)));
+    TermD = (Pp(TopRowWithSource(L_Q)));
+
+    L_TmLsGq0 = (TermD - TermA + TermB);
+    L_TmLsTmp = (TermC - TermB + TermA);
+
+    TermA = (Pp(Q_4d(R_Q)));
+    TermB = (Pm(Q_4d(R_Q)));
+    TermC = (Pm(TopRowWithSource(R_Q)));
+    TermD = (Pp(TopRowWithSource(R_Q)));
+
+    R_TmLsGq0 = (TermD - TermA + TermB);
+    R_TmLsTmp = (TermC - TermB + TermA);
+  }
+
+  std::vector<PropagatorField> R_TmLsGq(Ls,UGrid);
+  std::vector<PropagatorField> L_TmLsGq(Ls,UGrid);
+  for(int s=0;s<Ls;s++){
+    R_TmLsGq[s] = (Pm((R_Q)[(s)]) + Pp((R_Q)[((s)-1+Ls)%Ls]));
+    L_TmLsGq[s] = (Pm((L_Q)[(s)]) + Pp((L_Q)[((s)-1+Ls)%Ls]));
+  }
+
+  Gamma gmu=Gamma(Gmu[mu]);
+
+  q_out = Zero();
+  PropagatorField tmp(UGrid); 
+  for(int s=0;s<Ls;s++){
+
+    int sp = (s+1)%Ls;
+    int sr = Ls-1-s;
+    int srp= (sr+1)%Ls;
+
+    // Mobius parameters
+    auto b=this->bs[s];
+    auto c=this->cs[s];
+    auto bpc = 1.0/(b+c);  // -0.5 factor in gauge links
+    if (s == 0) {
+      p5d    =(b*Pm(L_TmLsGq[Ls-1])+ c*Pp(L_TmLsGq[Ls-1]) + b*Pp(L_TmLsTmp)   + c*Pm(L_TmLsTmp     ));
+      tmp    =(b*Pm(R_TmLsGq0)     + c*Pp(R_TmLsGq0     ) + b*Pp(R_TmLsGq[1]) + c*Pm(R_TmLsGq[1]));
+    } else if (s == Ls-1) {
+      p5d    =(b*Pm(L_TmLsGq0)     + c*Pp(L_TmLsGq0     ) + b*Pp(L_TmLsGq[1]) + c*Pm(L_TmLsGq[1]));
+      tmp    =(b*Pm(R_TmLsGq[Ls-1])+ c*Pp(R_TmLsGq[Ls-1]) + b*Pp(R_TmLsTmp)   + c*Pm(R_TmLsTmp   ));
+    } else {
+      p5d    =(b*Pm(L_TmLsGq[sr]) + c*Pp(L_TmLsGq[sr])+ b*Pp(L_TmLsGq[srp])+ c*Pm(L_TmLsGq[srp]));
+      tmp    =(b*Pm(R_TmLsGq[s])  + c*Pp(R_TmLsGq[s]) + b*Pp(R_TmLsGq[sp ])+ c*Pm(R_TmLsGq[sp]));
+    }
+    tmp    = Cshift(tmp,mu,1);
+    Impl::multLinkField(us_p5d,this->Umu,tmp,mu);
+    
+    gp5d=g5*p5d*g5;
+    gus_p5d=gmu*us_p5d;
+
+    C = bpc*(adj(gp5d)*us_p5d);
+    C-= bpc*(adj(gp5d)*gus_p5d);
+
+    if (s == 0) {
+      p5d    =(b*Pm(R_TmLsGq0)     + c*Pp(R_TmLsGq0  )    + b*Pp(R_TmLsGq[1]) + c*Pm(R_TmLsGq[1]));
+      tmp    =(b*Pm(L_TmLsGq[Ls-1])+ c*Pp(L_TmLsGq[Ls-1]) + b*Pp(L_TmLsTmp)   + c*Pm(L_TmLsTmp  ));
+    } else if (s == Ls-1) {
+      p5d    =(b*Pm(R_TmLsGq[Ls-1])+ c*Pp(R_TmLsGq[Ls-1]) + b*Pp(R_TmLsTmp)   + c*Pm(R_TmLsTmp  ));
+      tmp    =(b*Pm(L_TmLsGq0)     + c*Pp(L_TmLsGq0  )    + b*Pp(L_TmLsGq[1]) + c*Pm(L_TmLsGq[1]));
+    } else {
+      p5d    =(b*Pm(R_TmLsGq[s])  + c*Pp(R_TmLsGq[s])  + b*Pp(R_TmLsGq[sp ])+ c*Pm(R_TmLsGq[sp]));
+      tmp    =(b*Pm(L_TmLsGq[sr]) + c*Pp(L_TmLsGq[sr]) + b*Pp(L_TmLsGq[srp])+ c*Pm(L_TmLsGq[srp]));
+    }
+    tmp    = Cshift(tmp,mu,1);
+    Impl::multLinkField(us_p5d,this->Umu,tmp,mu);
+
+    gp5d=gmu*p5d;
+    gus_p5d=g5*us_p5d*g5;
+
+    C-= bpc*(adj(gus_p5d)*gp5d);
+    C-= bpc*(adj(gus_p5d)*p5d);
+
+    if (s < Ls/2) q_out += sgn*C;
+    else          q_out +=     C;
+    
+  }
+#endif
+}
+
+template <class Impl>
+void CayleyFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
+                                                PropagatorField &q_out,
+                                                PropagatorField &phys_src,
+                                                Current curr_type, 
+                                                unsigned int mu,
+                                                unsigned int tmin, 
+                                                unsigned int tmax,
+						ComplexField &ph)// Complex phase factor
+{
+  assert(mu>=0);
+  assert(mu<Nd);
+
+  int tshift = (mu == Nd-1) ? 1 : 0;
+
+#if 0
+  ////////////////////////////////////////////////
+  // SHAMIR CASE 
+  ////////////////////////////////////////////////
+  int Ls = this->Ls;
+  auto UGrid= this->GaugeGrid();
+  auto FGrid= this->FermionGrid();
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+  Gamma gmu=Gamma(Gmu[mu]);
+
+  PropagatorField L_Q(UGrid); 
+  PropagatorField R_Q(UGrid); 
+
+  PropagatorField tmp(UGrid);
+  PropagatorField Utmp(UGrid);
+  LatticeInteger zz (UGrid);   zz=0.0;
+  LatticeInteger lcoor(UGrid); LatticeCoordinate(lcoor,Nd-1);
+  for (int s=0;s<Ls;s++) {
+
+    RealD G_s = (curr_type == Current::Axial  ) ? ((s < Ls/2) ? -1 : 1) : 1;
+
+    ExtractSlice(R_Q, q_in, s , 0);
+
+    tmp    = Cshift(R_Q,mu,1);
+    Impl::multLinkField(Utmp,this->Umu,tmp,mu);
+    tmp    = G_s*( Utmp*ph - gmu*Utmp*ph ); // Forward hop
+    tmp    = where((lcoor>=tmin),tmp,zz); // Mask the time
+    tmp    = where((lcoor<=tmax),tmp,zz);
+    L_Q = tmp;
+
+    tmp    = R_Q*ph;
+    tmp    = Cshift(tmp,mu,-1);
+    Impl::multLinkField(Utmp,this->Umu,tmp,mu+Nd);// Adjoint link
+    tmp    = -G_s*( Utmp + gmu*Utmp ); 
+    tmp    = where((lcoor>=tmin+tshift),tmp,zz); // Mask the time 
+    tmp    = where((lcoor<=tmax+tshift),tmp,zz); // Position of current complicated
+    L_Q= L_Q+tmp;
+
+    InsertSlice(L_Q, q_out, s , 0);
+  }
+#endif
+
+#ifndef GRID_NVCC
+  ////////////////////////////////////////////////
+  // GENERAL CAYLEY CASE
+  ////////////////////////////////////////////////
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT,
+    Gamma::Algebra::Gamma5
+  };
+  Gamma gmu=Gamma(Gmu[mu]);
+  Gamma g5(Gamma::Algebra::Gamma5);
+
+  int Ls = this->Ls;
+  auto UGrid= this->GaugeGrid();
+  auto FGrid= this->FermionGrid();
+
+  std::vector<PropagatorField> R_Q(Ls,UGrid); 
+  PropagatorField L_Q(UGrid); 
+  PropagatorField tmp(UGrid);
+  PropagatorField Utmp(UGrid);
+
+  LatticeInteger zz (UGrid);   zz=0.0;
+  LatticeInteger lcoor(UGrid); LatticeCoordinate(lcoor,Nd-1);
+
+  for(int s=0;s<Ls;s++){
+    ExtractSlice(R_Q[s], q_in, s , 0);
+  }
+
+  PropagatorField R_TmLsGq0(UGrid); 
+  PropagatorField R_TmLsTmp(UGrid);
+  {
+    PropagatorField TermA(UGrid);
+    PropagatorField TermB(UGrid);
+    PropagatorField TermC(UGrid);
+    PropagatorField TermD(UGrid);
+
+    TermA = (Pp(Q_4d(R_Q)));
+    TermB = (Pm(Q_4d(R_Q)));
+    TermC = (Pm(TopRowWithSource(R_Q)));
+    TermD = (Pp(TopRowWithSource(R_Q)));
+
+    R_TmLsGq0 = (TermD - TermA + TermB);
+    R_TmLsTmp = (TermC - TermB + TermA);
+  }
+
+  std::vector<PropagatorField> R_TmLsGq(Ls,UGrid);
+  for(int s=0;s<Ls;s++){
+    R_TmLsGq[s] = (Pm((R_Q)[(s)]) + Pp((R_Q)[((s)-1+Ls)%Ls]));
+  }
+
+  std::vector<RealD> G_s(Ls,1.0);
+  if ( curr_type == Current::Axial ) {
+    for(int s=0;s<Ls/2;s++){
+      G_s[s] = -1.0;
+    }
+  }
+
+  for(int s=0;s<Ls;s++){
+
+    int sp = (s+1)%Ls;
+    int sr = Ls-1-s;
+    int srp= (sr+1)%Ls;
+
+    // Mobius parameters
+    auto b=this->bs[s];
+    auto c=this->cs[s];
+    //    auto bpc = G_s[s]*1.0/(b+c);  // -0.5 factor in gauge links
+
+    if (s == 0) {
+      tmp    =(b*Pm(R_TmLsGq0)     + c*Pp(R_TmLsGq0     ) + b*Pp(R_TmLsGq[1]) + c*Pm(R_TmLsGq[1]));
+    } else if (s == Ls-1) {
+      tmp    =(b*Pm(R_TmLsGq[Ls-1])+ c*Pp(R_TmLsGq[Ls-1]) + b*Pp(R_TmLsTmp)   + c*Pm(R_TmLsTmp   ));
+    } else {
+      tmp    =(b*Pm(R_TmLsGq[s])  + c*Pp(R_TmLsGq[s]) + b*Pp(R_TmLsGq[sp ])+ c*Pm(R_TmLsGq[sp]));
+    }
+
+    tmp    = Cshift(tmp,mu,1);
+    Impl::multLinkField(Utmp,this->Umu,tmp,mu);
+    tmp    = G_s[s]*( Utmp*ph - gmu*Utmp*ph ); // Forward hop
+    tmp    = where((lcoor>=tmin),tmp,zz); // Mask the time 
+    L_Q    = where((lcoor<=tmax),tmp,zz); // Position of current complicated
+
+    if (s == 0) {
+      tmp    =(b*Pm(R_TmLsGq0)     + c*Pp(R_TmLsGq0  )    + b*Pp(R_TmLsGq[1]) + c*Pm(R_TmLsGq[1]));
+    } else if (s == Ls-1) {
+      tmp    =(b*Pm(R_TmLsGq[Ls-1])+ c*Pp(R_TmLsGq[Ls-1]) + b*Pp(R_TmLsTmp)   + c*Pm(R_TmLsTmp  ));
+    } else {
+      tmp    =(b*Pm(R_TmLsGq[s])   + c*Pp(R_TmLsGq[s])    + b*Pp(R_TmLsGq[sp])+ c*Pm(R_TmLsGq[sp]));
+    }
+    tmp    = tmp *ph;
+    tmp    = Cshift(tmp,mu,-1);
+    Impl::multLinkField(Utmp,this->Umu,tmp,mu+Nd); // Adjoint link
+    tmp = -G_s[s]*( Utmp + gmu*Utmp );
+    tmp    = where((lcoor>=tmin+tshift),tmp,zz); // Mask the time 
+    L_Q   += where((lcoor<=tmax+tshift),tmp,zz); // Position of current complicated
+
+    InsertSlice(L_Q, q_out, s , 0);
+  }
+#endif
+}
+#undef Pp
+#undef Pm
+#undef Q_4d
+#undef TopRowWithSource
+
+
+
 #if 0
 template<class Impl>
 void CayleyFermion5D<Impl>::MooeeInternalCompute(int dag, int inv,

Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermion5DImplementation.h

@@ -611,6 +611,7 @@ template <class Impl>
 void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 								PropagatorField &q_in_2,
 								PropagatorField &q_out,
+								PropagatorField &src,
 								Current curr_type,
 								unsigned int mu)
 {
@@ -620,11 +621,12 @@ void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField
 template <class Impl>
 void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 							   PropagatorField &q_out,
+							   PropagatorField &src,
 							   Current curr_type,
 							   unsigned int mu, 
 							   unsigned int tmin,
-                                              unsigned int tmax,
+							   unsigned int tmax,
-					      ComplexField &lattice_cmplx)
+							   ComplexField &lattice_cmplx)
 {
   assert(0);
 

Grid/qcd/action/fermion/implementation/ImprovedStaggeredFermionImplementation.h

@@ -600,6 +600,7 @@ template <class Impl>
 void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 							      PropagatorField &q_in_2,
 							      PropagatorField &q_out,
+							      PropagatorField &src,
 							      Current curr_type,
 							      unsigned int mu)
 {
@@ -609,6 +610,7 @@ void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q
 template <class Impl>
 void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                                          PropagatorField &q_out,
+                                                         PropagatorField &src,
                                                          Current curr_type,
                                                          unsigned int mu, 
                                                          unsigned int tmin,

Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h

@@ -861,7 +861,6 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
  * Conserved current utilities for Wilson fermions, for contracting propagators
  * to make a conserved current sink or inserting the conserved current 
  * sequentially.
- ******************************************************************************/
 
 // Helper macro to reverse Simd vector. Fixme: slow, generic implementation.
 #define REVERSE_LS(qSite, qSiteRev, Nsimd) \
@@ -877,220 +876,10 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
     merge(qSiteRev, qSiteVec); \
 }
 
-//          psi = chiralProjectPlus(Result_s[Ls/2-1]);
+ ******************************************************************************/
-//          psi+= chiralProjectMinus(Result_s[Ls/2]);
-//         PJ5q+=localInnerProduct(psi,psi);
-
-template<class vobj> 
-Lattice<vobj> spProj5p(const Lattice<vobj> & in)
-{
-  GridBase *grid=in.Grid();
-  Gamma G5(Gamma::Algebra::Gamma5);
-  Lattice<vobj> ret(grid);
-  auto ret_v = ret.View();
-  auto in_v  =  in.View();
-  thread_for(ss,grid->oSites(),{
-    ret_v[ss] = in_v[ss] + G5*in_v[ss];
-  });
-  return ret;
-}
-template<class vobj> 
-Lattice<vobj> spProj5m(const Lattice<vobj> & in)
-{
-  Gamma G5(Gamma::Algebra::Gamma5);
-  GridBase *grid=in.Grid();
-  Lattice<vobj> ret(grid);
-  auto ret_v = ret.View();
-  auto in_v  =  in.View();
-  thread_for(ss,grid->oSites(),{
-    ret_v[ss] = in_v[ss] - G5*in_v[ss];
-  });
-  return ret;
-}
-
-template <class Impl>
-void WilsonFermion5D<Impl>::ContractJ5q(FermionField &q_in,ComplexField &J5q)
-{
-  conformable(GaugeGrid(), J5q.Grid());
-  conformable(q_in.Grid(), FermionGrid());
-
-  // 4d field
-  int Ls = this->Ls;
-  FermionField psi(GaugeGrid());
-  FermionField p_plus (GaugeGrid());
-  FermionField p_minus(GaugeGrid());
-  FermionField p(GaugeGrid());
-
-  ExtractSlice(p_plus , q_in, Ls/2   , 0);
-  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
-  p_plus = spProj5p(p_plus );
-  p_minus= spProj5m(p_minus);
-  p=p_plus+p_minus;
-  J5q = localInnerProduct(p,p);
-}
-
-template <class Impl>
-void WilsonFermion5D<Impl>::ContractJ5q(PropagatorField &q_in,ComplexField &J5q)
-{
-  conformable(GaugeGrid(), J5q.Grid());
-  conformable(q_in.Grid(), FermionGrid());
-
-  // 4d field
-  int Ls = this->Ls;
-  PropagatorField psi(GaugeGrid());
-  PropagatorField p_plus (GaugeGrid());
-  PropagatorField p_minus(GaugeGrid());
-  PropagatorField p(GaugeGrid());
-
-  ExtractSlice(p_plus , q_in, Ls/2   , 0);
-  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
-  p_plus = spProj5p(p_plus );
-  p_minus= spProj5m(p_minus);
-  p=p_plus+p_minus;
-  J5q = localInnerProduct(p,p);
-}
-
-template <class Impl>
-void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
-                                                     PropagatorField &q_in_2,
-                                                     PropagatorField &q_out,
-                                                     Current curr_type,
-                                                     unsigned int mu)
-{
-    conformable(q_in_1.Grid(), FermionGrid());
-    conformable(q_in_1.Grid(), q_in_2.Grid());
-    conformable(_FourDimGrid, q_out.Grid());
-
-    PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid());
-    unsigned int LLs = q_in_1.Grid()->_rdimensions[0];
-    q_out = Zero();
-
-    // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s), 
-    // q2(x + mu, Ls - 1 - s). 5D lattice so shift 4D coordinate mu by one.
-    tmp1 = Cshift(q_in_1, mu + 1, 1);
-    tmp2 = Cshift(q_in_2, mu + 1, 1);
-    auto q_in_1_v = q_in_1.View();
-    auto q_in_2_v = q_in_2.View();
-    auto tmp1_v   = tmp1.View();
-    auto tmp2_v   = tmp2.View();
-    auto q_out_v  = q_out.View();
-    auto Umu_v    = Umu.View();
-    thread_for(sU, Umu.Grid()->oSites(),{
-
-        unsigned int sF1 = sU * LLs;
-        unsigned int sF2 = (sU + 1) * LLs - 1;
-
-        for (unsigned int s = 0; s < LLs; ++s)
-        {
-            bool axial_sign = ((curr_type == Current::Axial) && \
-                               (s < (LLs / 2)));
-            SitePropagator qSite2, qmuSite2;
-
-            // If vectorised in 5th dimension, reverse q2 vector to match up
-            // sites correctly.
-            if (Impl::LsVectorised)
-            {
-                REVERSE_LS(q_in_2_v[sF2], qSite2, Ls / LLs);
-                REVERSE_LS(tmp2_v[sF2], qmuSite2, Ls / LLs);
-            }
-            else
-            {
-                qSite2   = q_in_2_v[sF2];
-                qmuSite2 = tmp2_v[sF2];
-            }
-            Kernels::ContractConservedCurrentSiteFwd(tmp1_v[sF1], 
-                                                     qSite2, 
-                                                     q_out_v[sU],
-                                                     Umu_v, sU, mu, axial_sign);
-            Kernels::ContractConservedCurrentSiteBwd(q_in_1_v[sF1],
-                                                     qmuSite2,
-                                                     q_out_v[sU],
-                                                     Umu_v, sU, mu, axial_sign);
-            sF1++;
-            sF2--;
-        }
-    });
-}
 
 
-template <class Impl>
-void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
-                                                PropagatorField &q_out,
-                                                Current curr_type, 
-                                                unsigned int mu,
-                                                unsigned int tmin, 
-                                                unsigned int tmax,
-						ComplexField &lattice_cmplx)
-{
-    conformable(q_in.Grid(), FermionGrid());
-    conformable(q_in.Grid(), q_out.Grid());
-    PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid());
-    unsigned int tshift = (mu == Tp) ? 1 : 0;
-    unsigned int LLs = q_in.Grid()->_rdimensions[0];
-    unsigned int LLt    = GridDefaultLatt()[Tp];
 
-    q_out = Zero();
-    LatticeInteger coords(_FourDimGrid);
-    LatticeCoordinate(coords, Tp);
-    
-    auto q_out_v = q_out.View();
-    auto tmp2_v  = tmp2.View();
-    auto coords_v= coords.View();
-    auto Umu_v   = Umu.View();
-    for (unsigned int s = 0; s < LLs; ++s)
-    {
-        bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
-	bool tadpole_sign = (curr_type == Current::Tadpole);
-	bool switch_sgn = tadpole_sign || axial_sign;
-
-
-        //forward direction: Need q(x + mu, s)*A(x)
-        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s) 
-        tmp = Cshift(tmp2, mu, 1);	 //q(x+mu,s)
-        tmp2 = tmp*lattice_cmplx;	 //q(x+mu,s)*A(x)	
-
-    	thread_for(sU, Umu.Grid()->oSites(),{
-            // Compute the sequential conserved current insertion only if our simd
-            // object contains a timeslice we need.
-            vPredicate t_mask;
-	    t_mask() = ((coords_v[sU] >= tmin) && (coords_v[sU] <= tmax));
-            Integer timeSlices = Reduce(t_mask());
-
-            if (timeSlices > 0)
-            {
-		unsigned int sF = sU * LLs + s;
-                Kernels::SeqConservedCurrentSiteFwd(tmp2_v[sU], 
-						    q_out_v[sF], Umu_v, sU,
-						    mu, t_mask, switch_sgn);
-            }
-
-        });
-
-        //backward direction: Need q(x - mu, s)*A(x-mu)
-        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s)
-        tmp = lattice_cmplx*tmp2;	 //q(x,s)*A(x)
-        tmp2 = Cshift(tmp, mu, -1);	 //q(x-mu,s)*A(x-mu,s)
-
-    	thread_for(sU, Umu.Grid()->oSites(),
-    	{
-	  vPredicate t_mask;
-	  t_mask()= ((coords_v[sU] >= (tmin + tshift)) && (coords_v[sU] <= (tmax + tshift)));
-
-	  //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3)	
-	  unsigned int t0 = 0;
-	  if((tmax==LLt-1) && (tshift==1)) t_mask() = (t_mask() || (coords_v[sU] == t0 ));
-	  
-	  Integer timeSlices = Reduce(t_mask());
-	  
-	  if (timeSlices > 0) {
-	    unsigned int sF = sU * LLs + s; 
-	    Kernels::SeqConservedCurrentSiteBwd(tmp2_v[sU], 
-						q_out_v[sF], Umu_v, sU,
-						mu, t_mask, axial_sign);
-	  }
-	});
-    }
-}
   
 NAMESPACE_END(Grid);
 

Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h

@@ -459,6 +459,7 @@ template <class Impl>
 void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
                                                    PropagatorField &q_in_2,
                                                    PropagatorField &q_out,
+                                                   PropagatorField &src,
                                                    Current curr_type,
                                                    unsigned int mu)
 {
@@ -466,6 +467,7 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
   conformable(_grid, q_in_1.Grid());
   conformable(_grid, q_in_2.Grid());
   conformable(_grid, q_out.Grid());
+#if 0
   PropagatorField tmp1(_grid), tmp2(_grid);
   q_out = Zero();
 
@@ -489,12 +491,15 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 					       q_out_v[sU],
 					       Umu_v, sU, mu);
   });
+#else
+#endif
 }
 
 
 template <class Impl>
 void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                               PropagatorField &q_out,
+                                              PropagatorField &src,
                                               Current curr_type,
                                               unsigned int mu,
                                               unsigned int tmin, 
@@ -503,6 +508,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 {
   conformable(_grid, q_in.Grid());
   conformable(_grid, q_out.Grid());
+#if 0
 
   //  Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
   Complex i(0.0,1.0);
@@ -556,6 +562,8 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 					  Umu_v, sU, mu, t_mask);
     }
   });
+#else
+#endif
 }
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h

@@ -493,131 +493,5 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
    assert(0 && " Kernel optimisation case not covered ");
   }
 
-/*******************************************************************************
- * Conserved current utilities for Wilson fermions, for contracting propagators
- * to make a conserved current sink or inserting the conserved current 
- * sequentially. Common to both 4D and 5D.
- ******************************************************************************/
-// N.B. Functions below assume a -1/2 factor within U.
-#define WilsonCurrentFwd(expr, mu) ((expr - Gamma::gmu[mu]*expr))
-#define WilsonCurrentBwd(expr, mu) ((expr + Gamma::gmu[mu]*expr))
-
-/*******************************************************************************
- * Name: ContractConservedCurrentSiteFwd
- * Operation: (1/2) * q2[x] * U(x) * (g[mu] - 1) * q1[x + mu]
- * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
- *        - Pass in q_in_1 shifted in +ve mu direction.
- ******************************************************************************/
-template<class Impl>
-void WilsonKernels<Impl>::ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1,
-							  const SitePropagator &q_in_2,
-							  SitePropagator &q_out,
-							  DoubledGaugeFieldView &U,
-							  unsigned int sU,
-							  unsigned int mu,
-							  bool switch_sign)
-{
-  SitePropagator result, tmp;
-  Gamma g5(Gamma::Algebra::Gamma5);
-
-  Impl::multLink(tmp, U[sU], q_in_1, mu);
-
-  result = g5 * adj(q_in_2) * g5 * WilsonCurrentFwd(tmp, mu);
-
-  if (switch_sign) {
-    q_out -= result;
-  } else {
-    q_out += result;
-  }
-}
-
-/*******************************************************************************
- * Name: ContractConservedCurrentSiteBwd
- * Operation: (1/2) * q2[x + mu] * U^dag(x) * (g[mu] + 1) * q1[x]
- * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
- *        - Pass in q_in_2 shifted in +ve mu direction.
- ******************************************************************************/
-template<class Impl>
-void WilsonKernels<Impl>::ContractConservedCurrentSiteBwd(const SitePropagator &q_in_1,
-							  const SitePropagator &q_in_2,
-							  SitePropagator &q_out,
-							  DoubledGaugeFieldView &U,
-							  unsigned int sU,
-							  unsigned int mu,
-							  bool switch_sign)
-{
-  SitePropagator result, tmp;
-  Gamma g5(Gamma::Algebra::Gamma5);
-
-  Impl::multLink(tmp, U[sU], q_in_1, mu + Nd);
-
-  result = g5 * adj(q_in_2) * g5 * WilsonCurrentBwd(tmp, mu);
-  if (switch_sign) {
-    q_out += result;
-  } else {
-    q_out -= result;
-  }
-}
-
-/*******************************************************************************
- * Name: SeqConservedCurrentSiteFwd
- * Operation: (1/2) * U(x) * (g[mu] - 1) * q[x + mu]
- * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
- *        - Pass in q_in shifted in +ve mu direction.
- ******************************************************************************/
-template<class Impl>
-void WilsonKernels<Impl>::SeqConservedCurrentSiteFwd(const SitePropagator &q_in,
-                                                     SitePropagator &q_out,
-                                                     DoubledGaugeFieldView &U,
-                                                     unsigned int sU,
-                                                     unsigned int mu,
-                                                     vPredicate t_mask,
-                                                     bool switch_sign)
-{
-  SitePropagator result;
-  
-  Impl::multLink(result, U[sU], q_in, mu);
-  result = WilsonCurrentFwd(result, mu);
-
-  // Zero any unwanted timeslice entries.
-  result = predicatedWhere(t_mask, result, 0.*result);
-  
-  if (switch_sign) {
-    q_out -= result;
-  } else {
-    q_out += result;
-  }
-}
-
-/*******************************************************************************
- * Name: SeqConservedCurrentSiteFwd
- * Operation: (1/2) * U^dag(x) * (g[mu] + 1) * q[x - mu]
- * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
- *        - Pass in q_in shifted in -ve mu direction.
- ******************************************************************************/
-template<class Impl>
-void WilsonKernels<Impl>::SeqConservedCurrentSiteBwd(const SitePropagator &q_in, 
-                                                     SitePropagator &q_out,
-                                                     DoubledGaugeFieldView &U,
-                                                     unsigned int sU,
-                                                     unsigned int mu,
-                                                     vPredicate t_mask,
-                                                     bool switch_sign)
-{
-  SitePropagator result;
-  Impl::multLink(result, U[sU], q_in, mu + Nd);
-  result = WilsonCurrentBwd(result, mu);
-
-  // Zero any unwanted timeslice entries.
-  result = predicatedWhere(t_mask, result, 0.*result);
-  
-  if (switch_sign) {
-    q_out += result;
-  } else {
-    q_out -= result;
-  }
-}
-
-
 NAMESPACE_END(Grid);
 

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/CayleyFermion5DInstantiationGparityWilsonImplD.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dcache.h>
-
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dvec.h>
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dgpu.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class CayleyFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/CayleyFermion5DInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/ContinuedFractionFermion5DInstantiationGparityWilsonImplD.cc

@@ -1,38 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/ContinuedFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class ContinuedFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/ContinuedFractionFermion5DInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/DomainWallEOFAFermionInstantiationGparityWilsonImplD.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class DomainWallEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/DomainWallEOFAFermionInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/MobiusEOFAFermionInstantiationGparityWilsonImplD.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/MobiusEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class MobiusEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/MobiusEOFAFermionInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/PartialFractionFermion5DInstantiationGparityWilsonImplD.cc

@@ -1,39 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/PartialFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/PartialFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
- 
-#include "impl.h"
-template class PartialFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/PartialFractionFermion5DInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonCloverFermionInstantiationGparityWilsonImplD.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
-
-    Copyright (C) 2017
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Guido Cossu <guido.cossu@ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/WilsonCloverFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonCloverFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonCloverFermionInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonFermion5DInstantiationGparityWilsonImplD.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonFermion5DInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonFermionInstantiationGparityWilsonImplD.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonFermionInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonKernelsInstantiationGparityWilsonImplD.cc

@@ -1,74 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandGparityImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-// Move these
-
-#include "impl.h"
-
-// G-parity requires more specialised implementation.
-template <>
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1, 
-							  const SitePropagator &q_in_2, 
-							  SitePropagator &q_out, 
-							  DoubledGaugeFieldView &U,
-							  unsigned int sU, 
-							  unsigned int mu, 
-							  bool switch_sign)
-{			
-  assert(0);		
-}			
-template <>		
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteBwd( const SitePropagator &q_in_1,
-							   const SitePropagator &q_in_2,
-							   SitePropagator &q_out,
-							   DoubledGaugeFieldView &U,
-							   unsigned int mu,
-							   unsigned int sU,
-							   bool switch_sign)
-{			
-  assert(0);		
-}
-
-HAND_SPECIALISE_GPARITY(IMPLEMENTATION);
-
-
-template class WilsonKernels<IMPLEMENTATION>; 
-
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonKernelsInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../WilsonKernelsInstantiationGparity.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonTMFermionInstantiationGparityWilsonImplD.cc

@@ -1,37 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/WilsonTMFermion.cc
-
-    Copyright (C) 2015
-
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/WilsonTMFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonTMFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonTMFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplD/WilsonTMFermionInstantiationGparityWilsonImplD.cc

@@ -0,0 +1 @@
+../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/CayleyFermion5DInstantiationGparityWilsonImplDF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dcache.h>
-
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dvec.h>
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dgpu.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class CayleyFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/CayleyFermion5DInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/ContinuedFractionFermion5DInstantiationGparityWilsonImplDF.cc

@@ -1,38 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/ContinuedFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class ContinuedFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/ContinuedFractionFermion5DInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/DomainWallEOFAFermionInstantiationGparityWilsonImplDF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class DomainWallEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/DomainWallEOFAFermionInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/MobiusEOFAFermionInstantiationGparityWilsonImplDF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/MobiusEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class MobiusEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/MobiusEOFAFermionInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/PartialFractionFermion5DInstantiationGparityWilsonImplDF.cc

@@ -1,39 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/PartialFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/PartialFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
- 
-#include "impl.h"
-template class PartialFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/PartialFractionFermion5DInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonCloverFermionInstantiationGparityWilsonImplDF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
-
-    Copyright (C) 2017
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Guido Cossu <guido.cossu@ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/WilsonCloverFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonCloverFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonCloverFermionInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonFermion5DInstantiationGparityWilsonImplDF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonFermion5DInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonFermionInstantiationGparityWilsonImplDF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonFermionInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonKernelsInstantiationGparityWilsonImplDF.cc

@@ -1,74 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandGparityImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-// Move these
-
-#include "impl.h"
-
-// G-parity requires more specialised implementation.
-template <>
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1, 
-							  const SitePropagator &q_in_2, 
-							  SitePropagator &q_out, 
-							  DoubledGaugeFieldView &U,
-							  unsigned int sU, 
-							  unsigned int mu, 
-							  bool switch_sign)
-{			
-  assert(0);		
-}			
-template <>		
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteBwd( const SitePropagator &q_in_1,
-							   const SitePropagator &q_in_2,
-							   SitePropagator &q_out,
-							   DoubledGaugeFieldView &U,
-							   unsigned int mu,
-							   unsigned int sU,
-							   bool switch_sign)
-{			
-  assert(0);		
-}
-
-HAND_SPECIALISE_GPARITY(IMPLEMENTATION);
-
-
-template class WilsonKernels<IMPLEMENTATION>; 
-
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonKernelsInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../WilsonKernelsInstantiationGparity.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonTMFermionInstantiationGparityWilsonImplDF.cc

@@ -1,37 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/WilsonTMFermion.cc
-
-    Copyright (C) 2015
-
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/WilsonTMFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonTMFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonTMFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplDF/WilsonTMFermionInstantiationGparityWilsonImplDF.cc

@@ -0,0 +1 @@
+../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/CayleyFermion5DInstantiationGparityWilsonImplF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dcache.h>
-
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dvec.h>
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dgpu.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class CayleyFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/CayleyFermion5DInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/ContinuedFractionFermion5DInstantiationGparityWilsonImplF.cc

@@ -1,38 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/ContinuedFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class ContinuedFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/ContinuedFractionFermion5DInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/DomainWallEOFAFermionInstantiationGparityWilsonImplF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class DomainWallEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/DomainWallEOFAFermionInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/MobiusEOFAFermionInstantiationGparityWilsonImplF.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/MobiusEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class MobiusEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/MobiusEOFAFermionInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../MobiusEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/PartialFractionFermion5DInstantiationGparityWilsonImplF.cc

@@ -1,39 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/PartialFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/PartialFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
- 
-#include "impl.h"
-template class PartialFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/PartialFractionFermion5DInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../PartialFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonCloverFermionInstantiationGparityWilsonImplF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/qcd/action/fermion/WilsonCloverFermion.cc
-
-    Copyright (C) 2017
-
-    Author: paboyle <paboyle@ph.ed.ac.uk>
-    Author: Guido Cossu <guido.cossu@ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-/*  END LEGAL */
-
-#include <Grid/Grid.h>
-#include <Grid/qcd/spin/Dirac.h>
-#include <Grid/qcd/action/fermion/WilsonCloverFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonCloverFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonCloverFermionInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonFermion5DInstantiationGparityWilsonImplF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonFermion5DInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../WilsonFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonFermionInstantiationGparityWilsonImplF.cc

@@ -1,40 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonFermionInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonKernelsInstantiationGparityWilsonImplF.cc

@@ -1,74 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsAsmImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonKernelsHandGparityImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-// Move these
-
-#include "impl.h"
-
-// G-parity requires more specialised implementation.
-template <>
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1, 
-							  const SitePropagator &q_in_2, 
-							  SitePropagator &q_out, 
-							  DoubledGaugeFieldView &U,
-							  unsigned int sU, 
-							  unsigned int mu, 
-							  bool switch_sign)
-{			
-  assert(0);		
-}			
-template <>		
-void WilsonKernels<IMPLEMENTATION>::ContractConservedCurrentSiteBwd( const SitePropagator &q_in_1,
-							   const SitePropagator &q_in_2,
-							   SitePropagator &q_out,
-							   DoubledGaugeFieldView &U,
-							   unsigned int mu,
-							   unsigned int sU,
-							   bool switch_sign)
-{			
-  assert(0);		
-}
-
-HAND_SPECIALISE_GPARITY(IMPLEMENTATION);
-
-
-template class WilsonKernels<IMPLEMENTATION>; 
-
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonKernelsInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../WilsonKernelsInstantiationGparity.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonTMFermionInstantiationGparityWilsonImplF.cc

@@ -1,37 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/WilsonTMFermion.cc
-
-    Copyright (C) 2015
-
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/WilsonTMFermion.h>
-#include <Grid/qcd/action/fermion/implementation/WilsonTMFermionImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class WilsonTMFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplF/WilsonTMFermionInstantiationGparityWilsonImplF.cc

@@ -0,0 +1 @@
+../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/CayleyFermion5DInstantiationGparityWilsonImplFH.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/WilsonKernels.cc
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dcache.h>
-
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dvec.h>
-			   //#include <Grid/qcd/action/fermion/implementation/CayleyFermion5Dgpu.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class CayleyFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);
-

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/CayleyFermion5DInstantiationGparityWilsonImplFH.cc

@@ -0,0 +1 @@
+../CayleyFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/ContinuedFractionFermion5DInstantiationGparityWilsonImplFH.cc

@@ -1,38 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/ContinuedFractionFermion5D.h>
-#include <Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class ContinuedFractionFermion5D<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/ContinuedFractionFermion5DInstantiationGparityWilsonImplFH.cc

@@ -0,0 +1 @@
+../ContinuedFractionFermion5DInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/DomainWallEOFAFermionInstantiationGparityWilsonImplFH.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/DomainWallEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class DomainWallEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/DomainWallEOFAFermionInstantiationGparityWilsonImplFH.cc

@@ -0,0 +1 @@
+../DomainWallEOFAFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/GparityWilsonImplFH/MobiusEOFAFermionInstantiationGparityWilsonImplFH.cc

@@ -1,44 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/fermion/MobiusEOFAFermion.cc
-
-Copyright (C) 2017
-
-Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-Author: David Murphy <dmurphy@phys.columbia.edu>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-			   /*  END LEGAL */
-
-#include <Grid/Grid_Eigen_Dense.h>
-#include <Grid/qcd/action/fermion/FermionCore.h>
-#include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionImplementation.h>
-#include <Grid/qcd/action/fermion/implementation/MobiusEOFAFermionCache.h>
-
-NAMESPACE_BEGIN(Grid);
-
-#include "impl.h"
-template class MobiusEOFAFermion<IMPLEMENTATION>; 
-
-NAMESPACE_END(Grid);