From e140b3f802549b63809fcf0ad88657f19497a75c Mon Sep 17 00:00:00 2001
From: David Murphy <djm2131@columbia.edu>
Date: Wed, 16 Aug 2017 23:36:23 -0400
Subject: [PATCH] Beginning to re-import Mobius EOFA

---
 lib/qcd/action/fermion/Fermion.h            |  31 +-
 lib/qcd/action/fermion/MobiusEOFAFermion.cc | 469 ++++++++++++++++++++
 lib/qcd/action/fermion/MobiusEOFAFermion.h  | 133 ++++++
 3 files changed, 631 insertions(+), 2 deletions(-)
 create mode 100644 lib/qcd/action/fermion/MobiusEOFAFermion.cc
 create mode 100644 lib/qcd/action/fermion/MobiusEOFAFermion.h

diff --git a/lib/qcd/action/fermion/Fermion.h b/lib/qcd/action/fermion/Fermion.h
index ac2a94b2..99d97895 100644
--- a/lib/qcd/action/fermion/Fermion.h
+++ b/lib/qcd/action/fermion/Fermion.h
@@ -1,6 +1,6 @@
 /*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+    Grid physics library, www.github.com/paboyle/Grid
 
     Source file: ./lib/qcd/action/fermion/Fermion_base_aggregate.h
 
@@ -38,6 +38,8 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 // - ContinuedFractionFermion5D.cc
 // - WilsonFermion.cc
 // - WilsonKernels.cc
+// - DomainWallEOFAFermion.cc
+// - MobiusEOFAFermion.cc
 //
 // The explicit instantiation is only avoidable if we move this source to headers and end up with include/parse/recompile
 // for EVERY .cc file. This define centralises the list and restores global push of impl cases
@@ -57,6 +59,7 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 #include <Grid/qcd/action/fermion/DomainWallFermion.h>
 #include <Grid/qcd/action/fermion/DomainWallEOFAFermion.h>
 #include <Grid/qcd/action/fermion/MobiusFermion.h>
+#include <Grid/qcd/action/fermion/MobiusEOFAFermion.h>
 #include <Grid/qcd/action/fermion/ZMobiusFermion.h>
 #include <Grid/qcd/action/fermion/SchurDiagTwoKappa.h>
 #include <Grid/qcd/action/fermion/ScaledShamirFermion.h>
@@ -129,6 +132,14 @@ typedef MobiusFermion<WilsonImplRL> MobiusFermionRL;
 typedef MobiusFermion<WilsonImplFH> MobiusFermionFH;
 typedef MobiusFermion<WilsonImplDF> MobiusFermionDF;
 
+typedef MobiusEOFAFermion<WilsonImplR> MobiusEOFAFermionR;
+typedef MobiusEOFAFermion<WilsonImplF> MobiusEOFAFermionF;
+typedef MobiusEOFAFermion<WilsonImplD> MobiusEOFAFermionD;
+
+typedef MobiusEOFAFermion<WilsonImplRL> MobiusEOFAFermionRL;
+typedef MobiusEOFAFermion<WilsonImplFH> MobiusEOFAFermionFH;
+typedef MobiusEOFAFermion<WilsonImplDF> MobiusEOFAFermionDF;
+
 typedef ZMobiusFermion<ZWilsonImplR> ZMobiusFermionR;
 typedef ZMobiusFermion<ZWilsonImplF> ZMobiusFermionF;
 typedef ZMobiusFermion<ZWilsonImplD> ZMobiusFermionD;
@@ -137,7 +148,7 @@ typedef ZMobiusFermion<ZWilsonImplRL> ZMobiusFermionRL;
 typedef ZMobiusFermion<ZWilsonImplFH> ZMobiusFermionFH;
 typedef ZMobiusFermion<ZWilsonImplDF> ZMobiusFermionDF;
 
-// Ls vectorised 
+// Ls vectorised
 typedef DomainWallFermion<DomainWallVec5dImplR> DomainWallFermionVec5dR;
 typedef DomainWallFermion<DomainWallVec5dImplF> DomainWallFermionVec5dF;
 typedef DomainWallFermion<DomainWallVec5dImplD> DomainWallFermionVec5dD;
@@ -162,6 +173,14 @@ typedef MobiusFermion<DomainWallVec5dImplRL> MobiusFermionVec5dRL;
 typedef MobiusFermion<DomainWallVec5dImplFH> MobiusFermionVec5dFH;
 typedef MobiusFermion<DomainWallVec5dImplDF> MobiusFermionVec5dDF;
 
+typedef MobiusEOFAFermion<DomainWallVec5dImplR> MobiusEOFAFermionVec5dR;
+typedef MobiusEOFAFermion<DomainWallVec5dImplF> MobiusEOFAFermionVec5dF;
+typedef MobiusEOFAFermion<DomainWallVec5dImplD> MobiusEOFAFermionVec5dD;
+
+typedef MobiusEOFAFermion<DomainWallVec5dImplRL> MobiusEOFAFermionVec5dRL;
+typedef MobiusEOFAFermion<DomainWallVec5dImplFH> MobiusEOFAFermionVec5dFH;
+typedef MobiusEOFAFermion<DomainWallVec5dImplDF> MobiusEOFAFermionVec5dDF;
+
 typedef ZMobiusFermion<ZDomainWallVec5dImplR> ZMobiusFermionVec5dR;
 typedef ZMobiusFermion<ZDomainWallVec5dImplF> ZMobiusFermionVec5dF;
 typedef ZMobiusFermion<ZDomainWallVec5dImplD> ZMobiusFermionVec5dD;
@@ -246,6 +265,14 @@ typedef MobiusFermion<GparityWilsonImplRL> GparityMobiusFermionRL;
 typedef MobiusFermion<GparityWilsonImplFH> GparityMobiusFermionFH;
 typedef MobiusFermion<GparityWilsonImplDF> GparityMobiusFermionDF;
 
+typedef MobiusEOFAFermion<GparityWilsonImplR> GparityMobiusFermionR;
+typedef MobiusEOFAFermion<GparityWilsonImplF> GparityMobiusFermionF;
+typedef MobiusEOFAFermion<GparityWilsonImplD> GparityMobiusFermionD;
+
+typedef MobiusEOFAFermion<GparityWilsonImplRL> GparityMobiusEOFAFermionRL;
+typedef MobiusEOFAFermion<GparityWilsonImplFH> GparityMobiusEOFAFermionFH;
+typedef MobiusEOFAFermion<GparityWilsonImplDF> GparityMobiusEOFAFermionDF;
+
 typedef ImprovedStaggeredFermion<StaggeredImplR> ImprovedStaggeredFermionR;
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;
 typedef ImprovedStaggeredFermion<StaggeredImplD> ImprovedStaggeredFermionD;
diff --git a/lib/qcd/action/fermion/MobiusEOFAFermion.cc b/lib/qcd/action/fermion/MobiusEOFAFermion.cc
new file mode 100644
index 00000000..9d0d8b11
--- /dev/null
+++ b/lib/qcd/action/fermion/MobiusEOFAFermion.cc
@@ -0,0 +1,469 @@
+/*************************************************************************************
+
+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>
+
+namespace Grid {
+namespace QCD {
+
+  template<class Impl>
+    MobiusEOFAFermion<Impl>::MobiusEOFAFermion(
+      GaugeField            &_Umu,
+      GridCartesian         &FiveDimGrid,
+      GridRedBlackCartesian &FiveDimRedBlackGrid,
+      GridCartesian         &FourDimGrid,
+      GridRedBlackCartesian &FourDimRedBlackGrid,
+      RealD _mq1, RealD _mq2, RealD _mq3,
+      RealD _shift, int _pm, RealD _M5,
+      RealD _b, RealD _c, const ImplParams &p) :
+    AbstractEOFAFermion<Impl>(_Umu, FiveDimGrid, FiveDimRedBlackGrid,
+        FourDimGrid, FourDimRedBlackGrid, _mq1, _mq2, _mq3,
+        _shift, _pm, _M5, _b, _c, p)
+    {
+      int Ls = this->Ls;
+
+      RealD eps = 1.0;
+      Approx::zolotarev_data *zdata = Approx::higham(eps, this->Ls);
+      assert(zdata->n == this->Ls);
+
+      std::cout << GridLogMessage << "MobiusEOFAFermion (b=" << _b <<
+        ",c=" << _c << ") with Ls=" << Ls << std::endl;
+      this->SetCoefficientsTanh(zdata, _b, _c);
+      std::cout << GridLogMessage << "EOFA parameters: (mq1=" << _mq1 <<
+        ",mq2=" << _mq2 << ",mq3=" << _mq3 << ",shift=" << _shift <<
+        ",pm=" << _pm << ")" << std::endl;
+
+      Approx::zolotarev_free(zdata);
+
+      if(_shift != 0.0){
+        SetCoefficientsPrecondShiftOps();
+      } else {
+        Mooee_shift.resize(Ls, 0.0);
+        MooeeInv_shift_lc.resize(Ls, 0.0);
+        MooeeInv_shift_norm.resize(Ls, 0.0);
+        MooeeInvDag_shift_lc.resize(Ls, 0.0);
+        MooeeInvDag_shift_norm.resize(Ls, 0.0);
+      }
+    }
+
+    /***************************************************************
+    /* Additional EOFA operators only called outside the inverter.
+    /* Since speed is not essential, simple axpby-style
+    /* implementations should be fine.
+    /***************************************************************/
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din, int sign, int dag)
+    {
+      int Ls = this->Ls;
+      RealD alpha = this->alpha;
+
+      Din = zero;
+      if((sign == 1) && (dag == 0)) { // \Omega_{+}
+        for(int s=0; s<Ls; ++s){
+          axpby_ssp(Din, 0.0, psi, 2.0*std::pow(1.0-alpha,Ls-s-1)/std::pow(1.0+alpha,Ls-s), psi, s, 0);
+        }
+      } else if((sign == -1) && (dag == 0)) { // \Omega_{-}
+        for(int s=0; s<Ls; ++s){
+          axpby_ssp(Din, 0.0, psi, 2.0*std::pow(1.0-alpha,s)/std::pow(1.0+alpha,s+1), psi, s, 0);
+        }
+      } else if((sign == 1 ) && (dag == 1)) { // \Omega_{+}^{\dagger}
+        for(int sp=0; sp<Ls; ++sp){
+          axpby_ssp(Din, 1.0, Din, 2.0*std::pow(1.0-alpha,Ls-sp-1)/std::pow(1.0+alpha,Ls-sp), psi, 0, sp);
+        }
+      } else if((sign == -1) && (dag == 1)) { // \Omega_{-}^{\dagger}
+        for(int sp=0; sp<Ls; ++sp){
+          axpby_ssp(Din, 1.0, Din, 2.0*std::pow(1.0-alpha,sp)/std::pow(1.0+alpha,sp+1), psi, 0, sp);
+        }
+      }
+    }
+
+    // This is the operator relating the usual Ddwf to TWQCD's EOFA Dirac operator (arXiv:1706.05843, Eqn. 6).
+    // It also relates the preconditioned and unpreconditioned systems described in Appendix B.2.
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::Dtilde(const FermionField& psi, FermionField& chi)
+    {
+      int Ls    = this->Ls;
+      RealD b   = 0.5 * ( 1.0 + this->alpha );
+      RealD c   = 0.5 * ( 1.0 - this->alpha );
+      RealD mq1 = this->mq1;
+
+      for(int s=0; s<Ls; ++s){
+        if(s == 0) {
+          axpby_ssp_pminus(chi, b, psi, -c, psi, s, s+1);
+          axpby_ssp_pplus (chi, 1.0, chi, mq1*c, psi, s, Ls-1);
+        } else if(s == (Ls-1)) {
+          axpby_ssp_pminus(chi, b, psi, mq1*c, psi, s, 0);
+          axpby_ssp_pplus (chi, 1.0, chi, -c, psi, s, s-1);
+        } else {
+          axpby_ssp_pminus(chi, b, psi, -c, psi, s, s+1);
+          axpby_ssp_pplus (chi, 1.0, chi, -c, psi, s, s-1);
+        }
+      }
+    }
+
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::DtildeInv(const FermionField& psi, FermionField& chi){ }
+
+    /*****************************************************************************************************/
+
+    template<class Impl>
+    RealD MobiusEOFAFermion<Impl>::M(const FermionField& psi, FermionField& chi)
+    {
+      int Ls = this->Ls;
+
+      FermionField Din(psi._grid);
+
+      this->Meooe5D(psi, Din);
+      this->DW(Din, chi, DaggerNo);
+      axpby(chi, 1.0, 1.0, chi, psi);
+      this->M5D(psi, chi);
+      return(norm2(chi));
+    }
+
+    template<class Impl>
+    RealD MobiusEOFAFermion<Impl>::Mdag(const FermionField& psi, FermionField& chi)
+    {
+      int Ls = this->Ls;
+
+      FermionField Din(psi._grid);
+
+      this->DW(psi, Din, DaggerYes);
+      this->MeooeDag5D(Din, chi);
+      this->M5Ddag(psi, chi);
+      axpby(chi, 1.0, 1.0, chi, psi);
+      return(norm2(chi));
+    }
+
+    /********************************************************************
+    /* Performance critical fermion operators called inside the inverter
+    /********************************************************************/
+
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::M5D(const FermionField& psi, FermionField& chi)
+    {
+        int Ls = this->Ls;
+
+        std::vector<Coeff_t> diag(Ls,1.0);
+        std::vector<Coeff_t> upper(Ls,-1.0);  upper[Ls-1] = this->mq1;
+        std::vector<Coeff_t> lower(Ls,-1.0);  lower[0]    = this->mq1;
+
+        // no shift term
+        if(this->shift == 0.0){ this->M5D(psi, chi, chi, lower, diag, upper); }
+
+        // fused M + shift operation
+        else{ this->M5D_shift(psi, chi, chi, lower, diag, upper, Mooee_shift); }
+    }
+
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField& psi, FermionField& chi)
+    {
+        int   Ls    = this->Ls;
+        int   pm    = this->pm;
+        RealD shift = this->shift;
+        RealD mq1   = this->mq1;
+        RealD mq2   = this->mq2;
+        RealD mq3   = this->mq3;
+
+        // coefficients for shift operator ( = shift*\gamma_{5}*R_{5}*\Delta_{\pm}(mq2,mq3)*P_{\pm} )
+        Coeff_t shiftp(0.0), shiftm(0.0);
+        if(shift != 0.0){
+          if(pm == 1){ shiftp = shift*(mq3-mq2); }
+          else{ shiftm = -shift*(mq3-mq2); }
+        }
+
+        std::vector<Coeff_t> diag(Ls,1.0);
+        std::vector<Coeff_t> upper(Ls,-1.0); upper[Ls-1] = mq1 + shiftp;
+        std::vector<Coeff_t> lower(Ls,-1.0); lower[0]    = mq1 + shiftm;
+
+        #if(0)
+            std::cout << GridLogMessage << "MobiusEOFAFermion::M5Ddag(FF&,FF&):" << std::endl;
+            for(int i=0; i<diag.size(); ++i){
+                std::cout << GridLogMessage << "diag[" << i << "] =" << diag[i] << std::endl;
+            }
+            for(int i=0; i<upper.size(); ++i){
+                std::cout << GridLogMessage << "upper[" << i << "] =" << upper[i] << std::endl;
+            }
+            for(int i=0; i<lower.size(); ++i){
+                std::cout << GridLogMessage << "lower[" << i << "] =" << lower[i] << std::endl;
+            }
+        #endif
+
+        this->M5Ddag(psi, chi, chi, lower, diag, upper);
+    }
+
+    // half checkerboard operations
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::Mooee(const FermionField& psi, FermionField& chi)
+    {
+        int Ls = this->Ls;
+
+        std::vector<Coeff_t> diag = this->bee;
+        std::vector<Coeff_t> upper(Ls);
+        std::vector<Coeff_t> lower(Ls);
+
+        for(int s=0; s<Ls; s++){
+          upper[s] = -this->cee[s];
+          lower[s] = -this->cee[s];
+        }
+        upper[Ls-1] = this->dm;
+        lower[0]    = this->dp;
+
+        this->M5D(psi, psi, chi, lower, diag, upper);
+    }
+
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::MooeeDag(const FermionField& psi, FermionField& chi)
+    {
+        int Ls = this->Ls;
+
+        std::vector<Coeff_t> diag = this->bee;
+        std::vector<Coeff_t> upper(Ls);
+        std::vector<Coeff_t> lower(Ls);
+
+        for(int s=0; s<Ls; s++){
+          upper[s] = -this->cee[s];
+          lower[s] = -this->cee[s];
+        }
+        upper[Ls-1] = this->dp;
+        lower[0]    = this->dm;
+
+        this->M5Ddag(psi, psi, chi, lower, diag, upper);
+    }
+
+    /****************************************************************************************/
+
+    //Zolo
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::SetCoefficientsInternal(RealD zolo_hi, std::vector<Coeff_t>& gamma, RealD b, RealD c)
+    {
+        int   Ls    = this->Ls;
+        int   pm    = this->pm;
+        RealD mq1   = this->mq1;
+        RealD mq2   = this->mq2;
+        RealD mq3   = this->mq3;
+        RealD shift = this->shift;
+
+        ////////////////////////////////////////////////////////
+        // Constants for the preconditioned matrix Cayley form
+        ////////////////////////////////////////////////////////
+        this->bs.resize(Ls);
+        this->cs.resize(Ls);
+        this->aee.resize(Ls);
+        this->aeo.resize(Ls);
+        this->bee.resize(Ls);
+        this->beo.resize(Ls);
+        this->cee.resize(Ls);
+        this->ceo.resize(Ls);
+
+        for(int i=0; i<Ls; ++i){
+          this->bee[i] = 4.0 - this->M5 + 1.0;
+          this->cee[i] = 1.0;
+        }
+
+        for(int i=0; i<Ls; ++i){
+          this->aee[i] = this->cee[i];
+          this->bs[i] = this->beo[i] = 1.0;
+          this->cs[i] = this->ceo[i] = 0.0;
+        }
+
+        //////////////////////////////////////////
+        // EOFA shift terms
+        //////////////////////////////////////////
+        if(pm == 1){
+          this->dp = mq1*this->cee[0] + shift*(mq3-mq2);
+          this->dm = mq1*this->cee[Ls-1];
+        } else if(this->pm == -1) {
+          this->dp = mq1*this->cee[0];
+          this->dm = mq1*this->cee[Ls-1] - shift*(mq3-mq2);
+        } else {
+          this->dp = mq1*this->cee[0];
+          this->dm = mq1*this->cee[Ls-1];
+        }
+
+        //////////////////////////////////////////
+        // LDU decomposition of eeoo
+        //////////////////////////////////////////
+        this->dee.resize(Ls+1);
+        this->lee.resize(Ls);
+        this->leem.resize(Ls);
+        this->uee.resize(Ls);
+        this->ueem.resize(Ls);
+
+        for(int i=0; i<Ls; ++i){
+
+          if(i < Ls-1){
+
+            this->lee[i] = -this->cee[i+1]/this->bee[i]; // sub-diag entry on the ith column
+
+            this->leem[i] = this->dm/this->bee[i];
+            for(int j=0; j<i; j++){ this->leem[i] *= this->aee[j]/this->bee[j]; }
+
+            this->dee[i] = this->bee[i];
+
+            this->uee[i] = -this->aee[i]/this->bee[i];   // up-diag entry on the ith row
+
+            this->ueem[i] = this->dp / this->bee[0];
+            for(int j=1; j<=i; j++){ this->ueem[i] *= this->cee[j]/this->bee[j]; }
+
+          } else {
+
+            this->lee[i]  = 0.0;
+            this->leem[i] = 0.0;
+            this->uee[i]  = 0.0;
+            this->ueem[i] = 0.0;
+
+          }
+        }
+
+        {
+          Coeff_t delta_d = 1.0 / this->bee[0];
+          for(int j=1; j<Ls-1; j++){ delta_d *= this->cee[j] / this->bee[j]; }
+          this->dee[Ls-1] = this->bee[Ls-1] + this->cee[0] * this->dm * delta_d;
+          this->dee[Ls] = this->bee[Ls-1] + this->cee[Ls-1] * this->dp * delta_d;
+        }
+
+        int inv = 1;
+        this->MooeeInternalCompute(0, inv, this->MatpInv, this->MatmInv);
+        this->MooeeInternalCompute(1, inv, this->MatpInvDag, this->MatmInvDag);
+    }
+
+    // Recompute Cayley-form coefficients for different shift
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::RefreshShiftCoefficients(RealD new_shift)
+    {
+        this->shift = new_shift;
+        Approx::zolotarev_data *zdata = Approx::higham(1.0, this->Ls);
+        this->SetCoefficientsTanh(zdata, 1.0, 0.0);
+    }
+
+    template<class Impl>
+    void MobiusEOFAFermion<Impl>::MooeeInternalCompute(int dag, int inv,
+        Vector<iSinglet<Simd> >& Matp, Vector<iSinglet<Simd> >& Matm)
+    {
+        int Ls = this->Ls;
+
+        GridBase* grid = this->FermionRedBlackGrid();
+        int LLs = grid->_rdimensions[0];
+
+        if(LLs == Ls){
+            return; // Not vectorised in 5th direction
+        }
+
+        Eigen::MatrixXcd Pplus  = Eigen::MatrixXcd::Zero(Ls,Ls);
+        Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls);
+
+        for(int s=0; s<Ls; s++){
+            Pplus(s,s)  = this->bee[s];
+            Pminus(s,s) = this->bee[s];
+        }
+
+        for(int s=0; s<Ls-1; s++){
+            Pminus(s,s+1) = -this->cee[s];
+        }
+
+        for(int s=0; s<Ls-1; s++){
+            Pplus(s+1,s) = -this->cee[s+1];
+        }
+
+        Pplus (0,Ls-1) = this->dp;
+        Pminus(Ls-1,0) = this->dm;
+
+        Eigen::MatrixXcd PplusMat ;
+        Eigen::MatrixXcd PminusMat;
+
+        #if(0)
+            std::cout << GridLogMessage << "Pplus:" << std::endl;
+            for(int s=0; s<Ls; ++s){
+                for(int ss=0; ss<Ls; ++ss){
+                    std::cout << Pplus(s,ss) << "\t";
+                }
+                std::cout << std::endl;
+            }
+            std::cout << GridLogMessage << "Pminus:" << std::endl;
+            for(int s=0; s<Ls; ++s){
+                for(int ss=0; ss<Ls; ++ss){
+                    std::cout << Pminus(s,ss) << "\t";
+                }
+                std::cout << std::endl;
+            }
+        #endif
+
+        if(inv) {
+            PplusMat  = Pplus.inverse();
+            PminusMat = Pminus.inverse();
+        } else {
+            PplusMat  = Pplus;
+            PminusMat = Pminus;
+        }
+
+        if(dag){
+            PplusMat.adjointInPlace();
+            PminusMat.adjointInPlace();
+        }
+
+        typedef typename SiteHalfSpinor::scalar_type scalar_type;
+        const int Nsimd = Simd::Nsimd();
+        Matp.resize(Ls*LLs);
+        Matm.resize(Ls*LLs);
+
+        for(int s2=0; s2<Ls; s2++){
+        for(int s1=0; s1<LLs; s1++){
+            int istride = LLs;
+            int ostride = 1;
+            Simd Vp;
+            Simd Vm;
+            scalar_type *sp = (scalar_type*) &Vp;
+            scalar_type *sm = (scalar_type*) &Vm;
+            for(int l=0; l<Nsimd; l++){
+                if(switcheroo<Coeff_t>::iscomplex()) {
+                    sp[l] = PplusMat (l*istride+s1*ostride,s2);
+                    sm[l] = PminusMat(l*istride+s1*ostride,s2);
+                } else {
+                    // if real
+                    scalar_type tmp;
+                    tmp = PplusMat (l*istride+s1*ostride,s2);
+                    sp[l] = scalar_type(tmp.real(),tmp.real());
+                    tmp = PminusMat(l*istride+s1*ostride,s2);
+                    sm[l] = scalar_type(tmp.real(),tmp.real());
+                }
+            }
+            Matp[LLs*s2+s1] = Vp;
+            Matm[LLs*s2+s1] = Vm;
+        }}
+    }
+
+    FermOpTemplateInstantiate(MobiusEOFAFermion);
+    GparityFermOpTemplateInstantiate(MobiusEOFAFermion);
+
+}}
diff --git a/lib/qcd/action/fermion/MobiusEOFAFermion.h b/lib/qcd/action/fermion/MobiusEOFAFermion.h
new file mode 100644
index 00000000..0a8d1788
--- /dev/null
+++ b/lib/qcd/action/fermion/MobiusEOFAFermion.h
@@ -0,0 +1,133 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/fermion/MobiusEOFAFermion.h
+
+Copyright (C) 2017
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: Peter Boyle <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 */
+#ifndef  GRID_QCD_MOBIUS_EOFA_FERMION_H
+#define  GRID_QCD_MOBIUS_EOFA_FERMION_H
+
+#include <Grid/qcd/action/fermion/AbstractEOFAFermion.h>
+
+namespace Grid {
+namespace QCD {
+
+  template<class Impl>
+  class MobiusEOFAFermion : public AbstractEOFAFermion<Impl>
+  {
+    public:
+      INHERIT_IMPL_TYPES(Impl);
+
+    public:
+      // Shift operator coefficients for red-black preconditioned Mobius EOFA
+      std::vector<Coeff_t> Mooee_shift;
+      std::vector<Coeff_t> MooeeInv_shift_lc;
+      std::vector<Coeff_t> MooeeInv_shift_norm;
+      std::vector<Coeff_t> MooeeInvDag_shift_lc;
+      std::vector<Coeff_t> MooeeInvDag_shift_norm;
+
+      virtual void Instantiatable(void) {};
+
+      // EOFA-specific operations
+      virtual void  Omega            (const FermionField& in, FermionField& out, int sign, int dag);
+      virtual void  Dtilde           (const FermionField& in, FermionField& out);
+      virtual void  DtildeInv        (const FermionField& in, FermionField& out);
+
+      // override multiply
+      virtual RealD M                (const FermionField& in, FermionField& out);
+      virtual RealD Mdag             (const FermionField& in, FermionField& out);
+
+      // half checkerboard operations
+      virtual void  Mooee            (const FermionField& in, FermionField& out);
+      virtual void  MooeeDag         (const FermionField& in, FermionField& out);
+      virtual void  MooeeInv         (const FermionField& in, FermionField& out);
+      virtual void  MooeeInv_shift   (const FermionField& in, FermionField& out);
+      virtual void  MooeeInvDag      (const FermionField& in, FermionField& out);
+      virtual void  MooeeInvDag_shift(const FermionField& in, FermionField& out);
+
+      virtual void   M5D             (const FermionField& psi, FermionField& chi);
+      virtual void   M5Ddag          (const FermionField& psi, FermionField& chi);
+
+      /////////////////////////////////////////////////////
+      // Instantiate different versions depending on Impl
+      /////////////////////////////////////////////////////
+      void M5D(const FermionField& psi, const FermionField& phi, FermionField& chi,
+        std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
+
+      void M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
+        std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
+        std::vector<Coeff_t>& shift_coeffs);
+
+      void M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi,
+        std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper);
+
+      void M5Ddag_shift(const FermionField& psi, const FermionField& phi, FermionField& chi,
+        std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper,
+        std::vector<Coeff_t>& shift_coeffs);
+
+      void MooeeInternal(const FermionField& in, FermionField& out, int dag, int inv);
+
+      void MooeeInternalCompute(int dag, int inv, Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
+
+      void MooeeInternalAsm(const FermionField& in, FermionField& out, int LLs, int site,
+        Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
+
+      void MooeeInternalZAsm(const FermionField& in, FermionField& out, int LLs, int site,
+        Vector<iSinglet<Simd>>& Matp, Vector<iSinglet<Simd>>& Matm);
+
+      virtual void RefreshShiftCoefficients(RealD new_shift);
+
+      // Constructors
+      MobiusEOFAFermion(GaugeField& _Umu, GridCartesian& FiveDimGrid, GridRedBlackCartesian& FiveDimRedBlackGrid,
+        GridCartesian& FourDimGrid, GridRedBlackCartesian& FourDimRedBlackGrid,
+        RealD _mq1, RealD _mq2, RealD _mq3, RealD _shift, int pm,
+        RealD _M5, RealD _b, RealD _c, const ImplParams& p=ImplParams());
+
+    protected:
+      void SetCoefficientsPrecondShiftOps(void);
+  };
+}}
+
+#define INSTANTIATE_DPERP_DWF_EOFA(A)\
+template void MobiusEOFAFermion<A>::M5D(const FermionField& psi, const FermionField& phi, FermionField& chi, \
+  std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
+template void MobiusEOFAFermion<A>::M5D_shift(const FermionField& psi, const FermionField& phi, FermionField& chi, \
+  std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper, std::vector<Coeff_t>& shift_coeffs); \
+template void MobiusEOFAFermion<A>::M5Ddag(const FermionField& psi, const FermionField& phi, FermionField& chi, \
+  std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper); \
+template void MobiusEOFAFermion<A>::M5Ddag_shift(const FermionField& psi, const FermionField& phi, FermionField& chi, \
+  std::vector<Coeff_t>& lower, std::vector<Coeff_t>& diag, std::vector<Coeff_t>& upper, std::vector<Coeff_t>& shift_coeffs); \
+template void MobiusEOFAFermion<A>::MooeeInv(const FermionField& psi, FermionField& chi); \
+template void MobiusEOFAFermion<A>::MooeeInv_shift(const FermionField& psi, FermionField& chi); \
+template void MobiusEOFAFermion<A>::MooeeInvDag(const FermionField& psi, FermionField& chi);
+template void MobiusEOFAFermion<A>::MooeeInvDag_shift(const FermionField& psi, FermionField& chi);
+
+#undef  MOBIUS_EOFA_DPERP_DENSE
+#define MOBIUS_EOFA_DPERP_CACHE
+#undef  MOBIUS_EOFA_DPERP_LINALG
+#define MOBIUS_EOFA_DPERP_VEC
+
+#endif