Reorganise a little to let the PV inverter be defined outside

the Reconstruct class.

This lets the multiple choices for PV inversion be composed without
changing the routine and no if/else case enumeration.

Implemented SchurDiagMooee PV inversion (red black) and Unprec PV inversion.
Red black cuts from 190 iterations to 90 iterations at 10^-12 on 8^4 test system

Will revisit multiple Schur options and add a Fourier based multishift PV inverse, similar
to the one Rudy Arthur did in BFM

Grid/algorithms/Algorithms.h

@@ -52,6 +52,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/CoarsenedMatrix.h>
 #include <Grid/algorithms/FFT.h>
 
+
 // EigCg
 // Pcg
 // Hdcg

Grid/algorithms/iterative/Reconstruct5Dprop.h

@@ -0,0 +1,186 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/SchurRedBlack.h
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+
+template<class Field>
+class PauliVillarsSolverUnprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverUnprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    MdagMLinearOperator<Matrix,Field> HermOp(_Matrix);
+
+    _Matrix.SetMass(1.0);
+    _Matrix.Mdag(src,A);
+    CG(HermOp,A,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field>
+class PauliVillarsSolverRBprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverRBprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    _Matrix.SetMass(1.0);
+    SchurRedBlackDiagMooeeSolve<Field> SchurSolver(CG);
+    SchurSolver(_Matrix,src,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class PVinverter> class Reconstruct5DfromPhysical {
+ private:
+  PVinverter & PauliVillarsSolver;
+ public:
+
+ /////////////////////////////////////////////////////
+ // First cut works, 10 Oct 2018.
+ //
+ // Must form a plan to get this into production for Zmobius acceleration
+ // of the Mobius exact AMA corrections.
+ //
+ // TODO : understand absence of contact term in eqns in Hantao's thesis
+ //        sol4 is contact term subtracted.
+ //
+ // Options
+ // a) Defect correction approach:
+ //    1) Compute defect from current soln (initially guess).
+ //       This is ...... outerToInner check !!!!
+ //    2) Deflated Zmobius solve to get 4d soln
+ //       Ensure deflation is working
+ //    3) Refine 5d Outer using the inner 4d delta soln
+ // 
+ // Step 1: localise PV inverse in a routine. [DONE]
+ // Step 2: Schur based PV inverse            [DONE]
+ // Step 3: Fourier accelerated PV inverse
+ // Step 4: 
+ /////////////////////////////////////////////////////
+ 
+  Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) 
+    : PauliVillarsSolver(_PauliVillarsSolver) 
+  { 
+  };
+
+
+   template<class Matrix>
+   void PV(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.M(src,sol);
+     _Matrix.SetMass(m);
+   }
+   template<class Matrix>
+   void PVdag(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.Mdag(src,sol);
+     _Matrix.SetMass(m);
+   }
+  template<class Matrix>
+  void operator() (Matrix & _Matrix,const Field &sol4,const Field &src4, Field &sol5){
+
+    int Ls =  _Matrix.Ls;
+
+    Field psi4(_Matrix.GaugeGrid());
+    Field psi(_Matrix.FermionGrid());
+    Field A  (_Matrix.FermionGrid());
+    Field B  (_Matrix.FermionGrid());
+    Field c  (_Matrix.FermionGrid());
+
+    typedef typename Matrix::Coeff_t Coeff_t;
+
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< " Reconstruct5Dprop: c.f. MADWF algorithm         " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    _Matrix.ImportPhysicalFermionSource(src4,B); 
+
+    ///////////////////////////////////////
+    // Set up c from src4
+    ///////////////////////////////////////
+    PauliVillarsSolver(_Matrix,B,A);
+    _Matrix.Pdag(A,c);
+
+    //////////////////////////////////////
+    // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+    //////////////////////////////////////
+    psi4 = - sol4;
+    InsertSlice(psi4, psi, 0   , 0);
+    for (int s=1;s<Ls;s++) {
+      ExtractSlice(psi4,c,s,0);
+       InsertSlice(psi4,psi,s,0);
+    }
+
+    /////////////////////////////
+    // Pdag PV^-1 Dm P 
+    /////////////////////////////
+    _Matrix.P(psi,B);
+    _Matrix.M(B,A);
+    PauliVillarsSolver(_Matrix,A,B);
+    _Matrix.Pdag(B,A);
+
+    //////////////////////////////
+    // Reinsert surface prop
+    //////////////////////////////
+    InsertSlice(sol4,A,0,0);
+
+    //////////////////////////////
+    // Convert from y back to x 
+    //////////////////////////////
+    _Matrix.P(A,sol5);
+    
+  }
+};
+
+}
+}

Grid/qcd/action/fermion/CayleyFermion5D.cc

@@ -68,6 +68,26 @@ void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &so
   ExtractSlice(exported4d, tmp, 0, 0);
 }
 template<class Impl>  
+void CayleyFermion5D<Impl>::P(const FermionField &psi, FermionField &chi)
+{
+  int Ls= this->Ls;
+  chi=zero;
+  for(int s=0;s<Ls;s++){
+    axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
+    axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s+1)%Ls);
+  }
+}
+template<class Impl>  
+void CayleyFermion5D<Impl>::Pdag(const FermionField &psi, FermionField &chi)
+{
+  int Ls= this->Ls;
+  chi=zero;
+  for(int s=0;s<Ls;s++){
+    axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
+    axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s-1+Ls)%Ls);
+  }
+}
+template<class Impl>  
 void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d)
 {
   int Ls = this->Ls;

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -93,6 +93,14 @@ namespace Grid {
       virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
       virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d);
 
+      ///////////////////////////////////////////////////////////////
+      // Support for MADWF tricks
+      ///////////////////////////////////////////////////////////////
+      RealD Mass(void) { return mass; };
+      void  SetMass(RealD _mass) { mass=_mass; } ;
+      void  P(const FermionField &psi, FermionField &chi);
+      void  Pdag(const FermionField &psi, FermionField &chi);
+
       /////////////////////////////////////////////////////
       // Instantiate different versions depending on Impl
       /////////////////////////////////////////////////////

Grid/qcd/action/fermion/FermionOperator.h

@@ -68,6 +68,7 @@ namespace Grid {
       virtual int    ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field
       virtual int    isTrivialEE(void) { return 0; };
       virtual RealD  Mass(void) {return 0.0;};
+      virtual void SetMass(RealD _mass) { return; };
 
       // half checkerboard operaions
       virtual void   Meooe       (const FermionField &in, FermionField &out)=0;

tests/debug/Test_cayley_cg.cc

@@ -27,6 +27,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/Reconstruct5Dprop.h>
 
 using namespace std;
 using namespace Grid;
@@ -75,6 +76,14 @@ void  TestCGprec(What & Ddwf,
 		 GridParallelRNG *RNG4,
 		 GridParallelRNG *RNG5);
 
+template<class What> 
+void  TestReconstruct5D(What & Ddwf, 
+			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			RealD mass, RealD M5,
+			GridParallelRNG *RNG4,
+			GridParallelRNG *RNG5);
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -152,6 +161,9 @@ void  TestCGinversions(What & Ddwf,
   TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
   std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl;
   TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
+
+  std::cout<<GridLogMessage << "Testing 5D PV reconstruction"<<std::endl;
+  TestReconstruct5D<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
 }
 
 template<class What> 
@@ -189,6 +201,53 @@ void  TestCGprec(What & Ddwf,
   CG(HermOpEO,src_o,result_o);
 }
 
+template<class What> 
+void  TestReconstruct5D(What & Ddwf, 
+			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			RealD mass, RealD M5,
+			GridParallelRNG *RNG4,
+			GridParallelRNG *RNG5)
+{
+  LatticeFermion src4   (UGrid); random(*RNG4,src4);
+  LatticeFermion res4   (UGrid); res4 = zero;
+
+  LatticeFermion src   (FGrid);
+  LatticeFermion src_NE(FGrid);
+  LatticeFermion result(FGrid);
+  LatticeFermion result_rec(FGrid);
+
+  MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
+  ConjugateGradient<LatticeFermion> CG(1.0e-12,10000);
+
+  Ddwf.ImportPhysicalFermionSource(src4,src);
+  Ddwf.Mdag(src,src_NE);
+  CG(HermOp,src_NE,result);
+
+  Ddwf.ExportPhysicalFermionSolution(result, res4);
+
+  Ddwf.M(result,src_NE);
+  src_NE = src_NE - src;
+  std::cout <<GridLogMessage<< " True residual is " << norm2(src_NE)<<std::endl;
+
+  std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
+
+  //  typedef PauliVillarsSolverUnprec<LatticeFermion> PVinverter;
+  typedef PauliVillarsSolverRBprec<LatticeFermion> PVinverter;
+  PVinverter PVinverse(CG);
+  Reconstruct5DfromPhysical<LatticeFermion,PVinverter> reconstructor(PVinverse);
+
+  reconstructor(Ddwf,res4,src4,result_rec);
+
+  std::cout <<GridLogMessage << "Result     "<<norm2(result)<<std::endl;
+  std::cout <<GridLogMessage << "Result_rec "<<norm2(result_rec)<<std::endl;
+
+  result_rec = result_rec - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
+
+}
+
+
 
 template<class What> 
 void  TestCGschur(What & Ddwf,