From 09288d633c9897916e20b85e70dccbf2b7a33733 Mon Sep 17 00:00:00 2001
From: Peter Boyle <paboyle@ph.ed.ac.uk>
Date: Sat, 10 Apr 2021 01:06:52 +0200
Subject: [PATCH] 4D pseudofermoin

---
 .../CayleyFermion5DImplementation.h           | 32 ++++++++++++++++++-
 1 file changed, 31 insertions(+), 1 deletion(-)

diff --git a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
index f11e9c44..8c01e064 100644
--- a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
+++ b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
@@ -112,7 +112,6 @@ void CayleyFermion5D<Impl>::ImportUnphysicalFermion(const FermionField &input4d,
   axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
   imported5d=tmp;
 }
-
 template<class Impl>  
 void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
 {
@@ -127,6 +126,37 @@ void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &inpu
   axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
   Dminus(tmp,imported5d);
 }
+////////////////////////////////////////////////////
+// Added for fourD pseudofermion det estimation
+////////////////////////////////////////////////////
+template<class Impl>  
+void CayleyFermion5D<Impl>::ImportFourDimPseudoFermion(const FermionField &input4d,FermionField &imported5d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  conformable(imported5d.Grid(),this->FermionGrid());
+  conformable(input4d.Grid()   ,this->GaugeGrid());
+  tmp = Zero();
+  InsertSlice(input4d, tmp, 0   , 0);
+  InsertSlice(input4d, tmp, Ls-1, 0);
+  axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, 0, 0);
+  axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
+  imported5d=tmp;
+}
+template<class Impl>  
+void CayleyFermion5D<Impl>::ExportFourDimPseudoFermion(const FermionField &solution5d,FermionField &exported4d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  tmp = solution5d;
+  conformable(solution5d.Grid(),this->FermionGrid());
+  conformable(exported4d.Grid(),this->GaugeGrid());
+  axpby_ssp_pminus(tmp, 0., solution5d, 1., solution5d, 0, 0);
+  axpby_ssp_pplus (tmp, 1., tmp       , 1., solution5d, 0, Ls-1);
+  ExtractSlice(exported4d, tmp, 0, 0);
+}
+
+// Dminus
 template<class Impl>  
 void CayleyFermion5D<Impl>::Dminus(const FermionField &psi, FermionField &chi)
 {