From f0229025e2e513f86258e5a25724446b4e7d3e3c Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Sat, 13 Oct 2018 19:55:03 +0100
Subject: [PATCH] MADWF working across a range of actions

---
 Grid/algorithms/iterative/ConjugateGradient.h |  14 +-
 Grid/qcd/action/fermion/CayleyFermion5D.cc    |   6 +-
 Grid/qcd/action/fermion/CayleyFermion5D.h     |  11 +-
 Grid/qcd/action/fermion/Fermion.h             |   1 +
 .../qcd/action/fermion/FourierAcceleratedPV.h |   2 +-
 Grid/qcd/action/fermion/MADWF.h               | 179 ++++++++++++++++++
 tests/debug/Test_cayley_cg.cc                 |  35 +++-
 7 files changed, 237 insertions(+), 11 deletions(-)
 create mode 100644 Grid/qcd/action/fermion/MADWF.h

diff --git a/Grid/algorithms/iterative/ConjugateGradient.h b/Grid/algorithms/iterative/ConjugateGradient.h
index ff4ba8ac..8a6a4bae 100644
--- a/Grid/algorithms/iterative/ConjugateGradient.h
+++ b/Grid/algorithms/iterative/ConjugateGradient.h
@@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
 	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;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
+	std::cout << GridLogPerformance << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 
diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.cc b/Grid/qcd/action/fermion/CayleyFermion5D.cc
index 0e0e7bee..a95ea4a0 100644
--- a/Grid/qcd/action/fermion/CayleyFermion5D.cc
+++ b/Grid/qcd/action/fermion/CayleyFermion5D.cc
@@ -485,9 +485,13 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
     
   double bpc = b+c;
   double bmc = b-c;
+  _b = b;
+  _c = c;
+  _gamma  = gamma; // Save the parameters so we can change mass later.
+  _zolo_hi= zolo_hi;
   for(int i=0; i < Ls; i++){
     as[i] = 1.0;
-    omega[i] = gamma[i]*zolo_hi; //NB reciprocal relative to Chroma NEF code
+    omega[i] = _gamma[i]*_zolo_hi; //NB reciprocal relative to Chroma NEF code
     assert(omega[i]!=Coeff_t(0.0));
     bs[i] = 0.5*(bpc/omega[i] + bmc);
     cs[i] = 0.5*(bpc/omega[i] - bmc);
diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h
index a84c0e8b..203a015e 100644
--- a/Grid/qcd/action/fermion/CayleyFermion5D.h
+++ b/Grid/qcd/action/fermion/CayleyFermion5D.h
@@ -97,7 +97,10 @@ namespace Grid {
       // Support for MADWF tricks
       ///////////////////////////////////////////////////////////////
       RealD Mass(void) { return mass; };
-      void  SetMass(RealD _mass) { mass=_mass; } ;
+      void  SetMass(RealD _mass) { 
+	mass=_mass; 
+	SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c);  // Reset coeffs
+      } ;
       void  P(const FermionField &psi, FermionField &chi);
       void  Pdag(const FermionField &psi, FermionField &chi);
 
@@ -147,6 +150,12 @@ namespace Grid {
       //    protected:
       RealD mass;
 
+      // Save arguments to SetCoefficientsInternal
+      std::vector<Coeff_t> _gamma;
+      RealD                _zolo_hi;
+      RealD                _b;
+      RealD                _c;
+
       // Cayley form Moebius (tanh and zolotarev)
       std::vector<Coeff_t> omega;
       std::vector<Coeff_t> bs;    // S dependent coeffs
diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h
index d0200ae7..77a4681f 100644
--- a/Grid/qcd/action/fermion/Fermion.h
+++ b/Grid/qcd/action/fermion/Fermion.h
@@ -91,6 +91,7 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 #include <Grid/qcd/action/fermion/FourierAcceleratedPV.h>
 #include <Grid/qcd/action/fermion/PauliVillarsInverters.h>
 #include <Grid/qcd/action/fermion/Reconstruct5Dprop.h>
+#include <Grid/qcd/action/fermion/MADWF.h>
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // More maintainable to maintain the following typedef list centrally, as more "impl" targets
diff --git a/Grid/qcd/action/fermion/FourierAcceleratedPV.h b/Grid/qcd/action/fermion/FourierAcceleratedPV.h
index d8441cf2..d6196eee 100644
--- a/Grid/qcd/action/fermion/FourierAcceleratedPV.h
+++ b/Grid/qcd/action/fermion/FourierAcceleratedPV.h
@@ -7,7 +7,7 @@
 
     Copyright (C) 2015
 
-Author: Christoph Lehner
+Author: Christoph Lehner (lifted with permission by Peter Boyle, brought back to Grid)
 Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 
     This program is free software; you can redistribute it and/or modify
diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h
new file mode 100644
index 00000000..f7ffea2f
--- /dev/null
+++ b/Grid/qcd/action/fermion/MADWF.h
@@ -0,0 +1,179 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/MADWF.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 Matrixo,class Matrixi,class PVinverter,class SchurSolver, class Guesser> 
+class MADWF 
+{
+ private:
+  typedef typename Matrixo::FermionField FermionFieldo;
+  typedef typename Matrixi::FermionField FermionFieldi;
+
+  PVinverter  & PauliVillarsSolvero;// For the outer field
+  SchurSolver & SchurSolveri;       // For the inner approx field
+  Guesser     & Guesseri;           // To deflate the inner approx solves
+
+  Matrixo & Mato;                   // Action object for outer
+  Matrixi & Mati;                   // Action object for inner
+
+  RealD target_resid;
+  int   maxiter;
+ public:
+
+  MADWF(Matrixo &_Mato,
+	Matrixi &_Mati, 
+	PVinverter &_PauliVillarsSolvero, 
+	SchurSolver &_SchurSolveri,
+	Guesser & _Guesseri,
+	RealD resid,
+	int _maxiter) :
+
+  Mato(_Mato),Mati(_Mati),
+    SchurSolveri(_SchurSolveri),
+    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri)
+  {   
+    target_resid=resid;
+    maxiter     =_maxiter; 
+  };
+
+  void operator() (const FermionFieldo &src4,FermionFieldo &sol5)
+  {
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< "  MADWF-like algorithm                           " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    FermionFieldo    c0(Mato.GaugeGrid()); // 4d 
+    FermionFieldo    y0(Mato.GaugeGrid()); // 4d
+
+    FermionFieldo    A(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    B(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    b(Mato.FermionGrid()); // 5d source
+
+    FermionFieldo    c(Mato.FermionGrid()); // PVinv source; reused so store
+    FermionFieldo    defect(Mato.FermionGrid()); // 5d source
+
+    FermionFieldi   ci(Mati.FermionGrid()); 
+    FermionFieldi   yi(Mati.FermionGrid()); 
+    FermionFieldi   xi(Mati.FermionGrid()); 
+    FermionFieldi srci(Mati.FermionGrid()); 
+    FermionFieldi   Ai(Mati.FermionGrid()); 
+
+    RealD m=Mati.Mass();
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    Mato.ImportPhysicalFermionSource(src4,b); 
+    std::cout << GridLogMessage << " src4 " <<norm2(src4)<<std::endl;
+    std::cout << GridLogMessage << " b    " <<norm2(b)<<std::endl;
+
+    defect = b;
+    sol5=zero;
+    for (int i=0;i<maxiter;i++) {
+
+      ///////////////////////////////////////
+      // Set up c0 from current defect
+      ///////////////////////////////////////
+      PauliVillarsSolvero(Mato,defect,A);
+      Mato.Pdag(A,c);
+      ExtractSlice(c0, c, 0 , 0);
+
+      ////////////////////////////////////////////////
+      // Solve the inner system with surface term c0
+      ////////////////////////////////////////////////
+      ci = zero;  
+      InsertSlice(c0,ci,0, 0);
+
+      // Dwm P y = Dwm x = D(1) P (c0,0,0,0)^T
+
+      Mati.P(ci,Ai);
+      Mati.SetMass(1.0);      Mati.M(Ai,srci);      Mati.SetMass(m);
+      SchurSolveri(Mati,srci,xi,Guesseri); 
+      Mati.Pdag(xi,yi);
+      ExtractSlice(y0, yi, 0 , 0);
+
+      //////////////////////////////////////
+      // Propagate solution back to outer system
+      // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+      //////////////////////////////////////
+      c0 = - y0;
+      InsertSlice(c0, c, 0   , 0);
+
+      /////////////////////////////
+      // Reconstruct the bulk solution Pdag PV^-1 Dm P 
+      /////////////////////////////
+      Mato.P(c,B);
+      Mato.M(B,A);
+      PauliVillarsSolvero(Mato,A,B);
+      Mato.Pdag(B,A);
+
+      //////////////////////////////
+      // Reinsert surface prop
+      //////////////////////////////
+      InsertSlice(y0,A,0,0);
+
+      //////////////////////////////
+      // Convert from y back to x 
+      //////////////////////////////
+      Mato.P(A,B);
+
+      //         sol5' = sol5 + M^-1 defect
+      //               = sol5 + M^-1 src - M^-1 M sol5  ...
+      sol5 = sol5 + B;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 update "<<std::endl;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 now "<<norm2(sol5)<<std::endl;
+      std::cout << GridLogMessage << " delta    "<<norm2(B)<<std::endl;
+
+       // New defect  = b - M sol5
+       Mato.M(sol5,A);
+       defect = b - A;
+
+       std::cout << GridLogMessage << " defect   "<<norm2(defect)<<std::endl;
+
+       double resid = ::sqrt(norm2(defect) / norm2(b));
+       std::cout << GridLogMessage << "Residual " << i << ": " << resid  << std::endl;
+       std::cout << GridLogMessage << "***************************************" <<std::endl;
+
+       if (resid < target_resid) {
+	 return;
+       }
+    }
+
+    std::cout << GridLogMessage << "MADWF : Exceeded maxiter "<<std::endl;
+    assert(0);
+
+  }
+
+};
+
+}}
diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc
index fbc978a6..361b07e9 100644
--- a/tests/debug/Test_cayley_cg.cc
+++ b/tests/debug/Test_cayley_cg.cc
@@ -250,10 +250,13 @@ void  TestReconstruct5D(What & Ddwf,
   LatticeFermion src_NE(FGrid);
   LatticeFermion result(FGrid);
   LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
 
   MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
   double Resid = 1.0e-12;
+  double Residi = 1.0e-6;
   ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermion> CGi(Residi,10000);
 
   Ddwf.ImportPhysicalFermionSource(src4,src);
   Ddwf.Mdag(src,src_NE);
@@ -271,7 +274,8 @@ void  TestReconstruct5D(What & Ddwf,
   // RBprec PV inverse
   ////////////////////////////
   typedef LatticeFermion Field;
-  typedef SchurRedBlackDiagMooeeSolve<Field> SchurSolverType; 
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverType; 
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverTypei; 
   typedef PauliVillarsSolverRBprec<Field,SchurSolverType> PVinverter;
   SchurSolverType SchurSolver(CG);
   PVinverter      PVinverse(SchurSolver);
@@ -286,6 +290,19 @@ void  TestReconstruct5D(What & Ddwf,
   result_rec = result_rec - result;
   std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
 
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolveri(CGi);
+  ZeroGuesser<LatticeFermion> Guess;
+  MADWF<What,What,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermion> > 
+    madwf(Ddwf,Ddwf,PVinverse,SchurSolveri,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
+
 }
 template<class What> 
 void  TestReconstruct5DFA(What & Ddwf, 
@@ -303,10 +320,13 @@ void  TestReconstruct5DFA(What & Ddwf,
   LatticeFermion src_NE(FGrid);
   LatticeFermion result(FGrid);
   LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
 
   MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
   double Resid = 1.0e-12;
+  double Residi = 1.0e-5;
   ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermion> CGi(Residi,10000);
 
   Ddwf.ImportPhysicalFermionSource(src4,src);
   Ddwf.Mdag(src,src_NE);
@@ -324,6 +344,7 @@ void  TestReconstruct5DFA(What & Ddwf,
   // Fourier accel PV inverse
   ////////////////////////////
   typedef LatticeFermion Field;
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverTypei; 
   typedef PauliVillarsSolverFourierAccel<LatticeFermion,LatticeGaugeField> PVinverter;
   PVinverter PVinverse(Umu,CG);
 
@@ -337,6 +358,18 @@ void  TestReconstruct5DFA(What & Ddwf,
   result_rec = result_rec - result;
   std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
 
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolver(CGi);
+  ZeroGuesser<LatticeFermion> Guess;
+  MADWF<What,What,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermion> > 
+    madwf(Ddwf,Ddwf,PVinverse,SchurSolver,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
 }