From 751fae9f0d64306944dc5f01b055b230a9a083e2 Mon Sep 17 00:00:00 2001
From: Chulwoo Jung <chulwoo@bnl.gov>
Date: Tue, 10 Jul 2018 12:18:12 -0700
Subject: [PATCH] Changing boundary phase to be always double

---
 extras/Hadrons/Modules/MAction/DWF.hpp        |  2 +-
 extras/Hadrons/Modules/MAction/Wilson.hpp     |  2 +-
 .../Hadrons/Modules/MAction/WilsonClover.hpp  |  2 +-
 lib/qcd/action/ActionParams.h                 |  4 +--
 tests/solver/Test_mobius_bcg.cc               | 27 ++++++++++++-------
 5 files changed, 23 insertions(+), 14 deletions(-)

diff --git a/extras/Hadrons/Modules/MAction/DWF.hpp b/extras/Hadrons/Modules/MAction/DWF.hpp
index 4dfd06cf..00b16d10 100644
--- a/extras/Hadrons/Modules/MAction/DWF.hpp
+++ b/extras/Hadrons/Modules/MAction/DWF.hpp
@@ -117,7 +117,7 @@ void TDWF<FImpl>::setup(void)
     auto &grb4 = *env().getRbGrid();
     auto &g5   = *env().getGrid(par().Ls);
     auto &grb5 = *env().getRbGrid(par().Ls);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    std::vector<ComplexD> boundary = strToVec<ComplexD>(par().boundary);
     typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
diff --git a/extras/Hadrons/Modules/MAction/Wilson.hpp b/extras/Hadrons/Modules/MAction/Wilson.hpp
index 6467b3ee..77cd0403 100644
--- a/extras/Hadrons/Modules/MAction/Wilson.hpp
+++ b/extras/Hadrons/Modules/MAction/Wilson.hpp
@@ -110,7 +110,7 @@ void TWilson<FImpl>::setup(void)
     auto &U      = envGet(LatticeGaugeField, par().gauge);
     auto &grid   = *env().getGrid();
     auto &gridRb = *env().getRbGrid();
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    std::vector<ComplexD> boundary = strToVec<ComplexD>(par().boundary);
     typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
diff --git a/extras/Hadrons/Modules/MAction/WilsonClover.hpp b/extras/Hadrons/Modules/MAction/WilsonClover.hpp
index c369f086..039a8b84 100644
--- a/extras/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/extras/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -121,7 +121,7 @@ void TWilsonClover<FImpl>::setup(void)
     auto &U      = envGet(LatticeGaugeField, par().gauge);
     auto &grid   = *env().getGrid();
     auto &gridRb = *env().getRbGrid();
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    std::vector<ComplexD> boundary = strToVec<ComplexD>(par().boundary);
     typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid, gridRb, par().mass,
 						  par().csw_r,
diff --git a/lib/qcd/action/ActionParams.h b/lib/qcd/action/ActionParams.h
index d25b60a9..8d4dcf49 100644
--- a/lib/qcd/action/ActionParams.h
+++ b/lib/qcd/action/ActionParams.h
@@ -44,11 +44,11 @@ namespace QCD {
   
   struct WilsonImplParams {
     bool overlapCommsCompute;
-    std::vector<Complex> boundary_phases;
+    std::vector<ComplexD> boundary_phases;
     WilsonImplParams() : overlapCommsCompute(false) {
       boundary_phases.resize(Nd, 1.0);
     };
-    WilsonImplParams(const std::vector<Complex> phi)
+    WilsonImplParams(const std::vector<ComplexD> phi)
       : boundary_phases(phi), overlapCommsCompute(false) {}
   };
 
diff --git a/tests/solver/Test_mobius_bcg.cc b/tests/solver/Test_mobius_bcg.cc
index 4b22fb18..f56c73ca 100644
--- a/tests/solver/Test_mobius_bcg.cc
+++ b/tests/solver/Test_mobius_bcg.cc
@@ -38,7 +38,7 @@ int main (int argc, char ** argv)
   typedef typename MobiusFermionR::ComplexField ComplexField; 
   typename MobiusFermionR::ImplParams params; 
 
-  const int Ls=24;
+  const int Ls=12;
 
   Grid_init(&argc,&argv);
 
@@ -49,6 +49,10 @@ int main (int argc, char ** argv)
   std::vector<int> split_coor (mpi_layout.size(),1);
   std::vector<int> split_dim (mpi_layout.size(),1);
 
+  std::vector<ComplexD> boundary_phases(Nd,1.);
+  boundary_phases[Nd-1]=-1.;
+  params.boundary_phases = boundary_phases;
+
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
 								   GridDefaultSimd(Nd,vComplex::Nsimd()),
 								   GridDefaultMpi());
@@ -139,7 +143,7 @@ int main (int argc, char ** argv)
 
   LatticeGaugeField Umu(UGrid); 
   FieldMetaData header;
-    std::string file("./lat.in.24ID");
+    std::string file("./lat.in.32IDfine");
     SU3::ColdConfiguration(Umu);
     std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
   if(1) { 
@@ -203,14 +207,14 @@ int main (int argc, char ** argv)
   // Set up N-solvers as trivially parallel
   ///////////////////////////////////////////////////////////////
   std::cout << GridLogMessage << " Building the solvers"<<std::endl;
-//  RealD mass=0.00107;
-  RealD mass=0.01;
+  RealD mass=0.00107;
+//  RealD mass=0.01;
   RealD M5=1.8;
-  RealD mobius_factor=4;
+  RealD mobius_factor=32./12.;
   RealD mobius_b=0.5*(mobius_factor+1.);
   RealD mobius_c=0.5*(mobius_factor-1.);
-  MobiusFermionR Dchk(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,mobius_b,mobius_c);
-  MobiusFermionR Ddwf(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5,mobius_b,mobius_c);
+  MobiusFermionR Dchk(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,mobius_b,mobius_c,params);
+  MobiusFermionR Ddwf(s_Umu,*SFGrid,*SFrbGrid,*SGrid,*SrbGrid,mass,M5,mobius_b,mobius_c,params);
 
   std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
   std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
@@ -218,8 +222,9 @@ int main (int argc, char ** argv)
 
   MdagMLinearOperator<MobiusFermionR,FermionField> HermOp(Ddwf);
   MdagMLinearOperator<MobiusFermionR,FermionField> HermOpCk(Dchk);
-  ConjugateGradient<FermionField> CG((stp),10000);
+  ConjugateGradient<FermionField> CG((stp),100000);
   s_res = zero;
+if(0){
 //  CG(HermOp,s_src,s_res);
 
   std::cout << GridLogMessage << " split residual norm "<<norm2(s_res)<<std::endl;
@@ -247,17 +252,21 @@ int main (int argc, char ** argv)
     HermOpCk.HermOp(result[n],tmp); tmp = tmp - src[n];
     std::cout << GridLogMessage<<" resid["<<n<<"]  "<< norm2(tmp)/norm2(src[n])<<std::endl;
   }
+}
 
 // faking enlarged/cooperative CG
+if(0){
+    std::cout << GridLogMessage<<" Trying blocking enlarged CG" <<std::endl;
   assert(me < nrhs);
   if (me>0) src[me] = src[0];
   for(int s=0;s<nrhs;s++){
      result[s]=zero;
      if(s!=me) src[s] = zero;
   }
+}
 
   int blockDim = 0;//not used for BlockCGVec
-  BlockConjugateGradient<FermionField>    BCGV  (BlockCGVec,blockDim,stp,10000);
+  BlockConjugateGradient<FermionField>    BCGV  (BlockCGVec,blockDim,stp,100000);
   BCGV.PrintInterval=10;
 {
   BCGV(HermOpCk,src,result);