Changing boundary phase to be always double

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);

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);

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,

lib/algorithms/Algorithms.h

@@ -49,7 +49,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/algorithms/iterative/BlockConjugateGradient.h>
 #include <Grid/algorithms/iterative/ConjugateGradientReliableUpdate.h>
 #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
-#include <Grid/algorithms/iterative/InexactPrecConjugateGradient.h>
 #include <Grid/algorithms/CoarsenedMatrix.h>
 #include <Grid/algorithms/FFT.h>
 

lib/algorithms/approx/MultiShiftFunction.cc

@@ -53,10 +53,4 @@ void MultiShiftFunction::csv(std::ostream &out)
   }
   return;
 }
-
-RealD __InverseApproximation(RealD x)
-{
-  return 1.0/x;
-}
-
 }

lib/algorithms/approx/MultiShiftFunction.h

@@ -63,7 +63,5 @@ public:
   }
 
 };
-RealD __InverseApproximation(RealD x);
-
 }
 #endif

lib/algorithms/iterative/BlockConjugateGradient.h

@@ -33,7 +33,7 @@ directory
 
 namespace Grid {
 
-enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS };
+enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec };
 
 //////////////////////////////////////////////////////////////////////////
 // Block conjugate gradient. Dimension zero should be the block direction
@@ -54,9 +54,10 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
   RealD Tolerance;
   Integer MaxIterations;
   Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
+  Integer PrintInterval; //GridLogMessages or Iterative
   
   BlockConjugateGradient(BlockCGtype cgtype,int _Orthog,RealD tol, Integer maxit, bool err_on_no_conv = true)
-    : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv)
+    : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv),PrintInterval(100)
   {};
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -127,6 +128,14 @@ void operator()(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
     assert(0);
   }
 }
+void operator()(LinearOperatorBase<Field> &Linop, const std::vector<Field> &Src, std::vector<Field> &Psi) 
+{
+  if ( CGtype == BlockCGVec ) {
+    BlockCGVecsolve(Linop,Src,Psi);
+  } else {
+    assert(0);
+  }
+}
 
 ////////////////////////////////////////////////////////////////////////////
 // BlockCGrQ implementation:
@@ -600,6 +609,272 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
   IterationsToComplete = k;
 }
 
+void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector<Field> &X, std::vector<Field> &Y){
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<Nblock;bp++) {
+    m(b,bp) = innerProduct(X[b],Y[bp]);  
+    }
+}
+double HermCheck( Eigen::MatrixXcd &m,  const std::string &str, int ForceHerm=1 , int Print = 0) {
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<=b;bp++) {
+    if(Print)
+    std::cout<<GridLogMessage << "HermCheck "<<str<<" "<<b<<" "<<bp<<" : "<< m(b,bp) <<" "<<conj(m(bp,b))<<" " <<m(b,bp)-conj(m(bp,b)) <<std::endl;
+    if(ForceHerm){
+      if(b==bp) m(b,b) = real(m(b,b));
+      else{
+        auto temp = 0.5*(m(b,bp)+conj(m(bp,b)));
+	m(b,bp) = temp;
+        m(bp,b) = conj(temp);
+      }
+    }
+  }
+}
+
+void BlockCGVecsolve(LinearOperatorBase<Field> &Linop, const std::vector<Field> &Src, std::vector<Field> &Psi) 
+{
+//  int Orthog = blockDim; // First dimension is block dim; this is an assumption
+//  Nblock = Src._grid->_fdimensions[Orthog];
+  Nblock = Src.size();
+  assert(Nblock == Psi.size());
+
+  std::cout<<GridLogMessage<<" Block Conjugate Gradient :  Nblock "<<Nblock<<std::endl;
+
+  for(int b=0;b<Nblock;b++){
+  Psi[b].checkerboard = Src[0].checkerboard;
+  conformable(Psi[b], Src[b]);
+  }
+
+  Field Fake(Src[0]);
+
+  std::vector<Field> P(Nblock,Fake);
+// P.resize(Nblock);
+  std::vector<Field> AP(Nblock,Fake); 
+//AP.resize(Nblock);
+  std::vector<Field> R(Nblock,Fake); 
+  std::vector<Field> TMP(Nblock,Fake); 
+//R.resize(Nblock);
+  
+  Eigen::MatrixXcd m_pAp    = Eigen::MatrixXcd::Identity(Nblock,Nblock);
+  Eigen::MatrixXcd m_pAp_inv= Eigen::MatrixXcd::Identity(Nblock,Nblock);
+  Eigen::MatrixXcd m_rr     = Eigen::MatrixXcd::Zero(Nblock,Nblock);
+  Eigen::MatrixXcd m_rr_inv = Eigen::MatrixXcd::Zero(Nblock,Nblock);
+
+  Eigen::MatrixXcd m_alpha      = Eigen::MatrixXcd::Zero(Nblock,Nblock);
+  Eigen::MatrixXcd m_beta   = Eigen::MatrixXcd::Zero(Nblock,Nblock);
+
+  // Initial residual computation & set up
+  std::vector<RealD> residuals(Nblock);
+  std::vector<RealD> ssq(Nblock);
+
+//  sliceNorm(ssq,Src,Orthog);
+  for(int b=0;b<Nblock;b++){ ssq[b] = norm2(Src[b]);}
+  RealD sssum=0;
+  for(int b=0;b<Nblock;b++) sssum+=ssq[b];
+
+//  sliceNorm(residuals,Src,Orthog);
+  for(int b=0;b<Nblock;b++){ residuals[b] = norm2(Src[b]);}
+  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
+
+//  sliceNorm(residuals,Psi,Orthog);
+  for(int b=0;b<Nblock;b++){ residuals[b] = norm2(Psi[b]);}
+  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
+
+  // Initial search dir is guess
+  for(int b=0;b<Nblock;b++) Linop.HermOp(Psi[b], AP[b]);
+  for(int b=0;b<Nblock;b++) 
+  std::cout << b << " Psi " << norm2(Psi[b]) <<" AP "<<norm2(AP[b])<<std::endl;
+  
+
+  /************************************************************************
+   * Block conjugate gradient (Stephen Pickles, thesis 1995, pp 71, O Leary 1980)
+   ************************************************************************
+   * O'Leary : R = B - A X
+   * O'Leary : P = M R ; preconditioner M = 1
+   * O'Leary : alpha = PAP^{-1} RMR
+   * O'Leary : beta  = RMR^{-1}_old RMR_new
+   * O'Leary : X=X+Palpha
+   * O'Leary : R_new=R_old-AP alpha
+   * O'Leary : P=MR_new+P beta
+   */
+
+  for(int b=0;b<Nblock;b++){
+  R[b] = Src[b] - AP[b]; //R_0  
+  P[b] = R[b];  // P_1
+  }
+//  sliceInnerProductMatrix(m_rr,R,R,Orthog);
+  InnerProductMatrix(m_rr,R,R);
+  HermCheck(m_rr, "R_0 R_0",1,1);
+  HermCheck(m_rr, "R_0 R_0",0,1);
+#if 0
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<Nblock;bp++) {
+    m_rr(b,bp) = innerProduct(R[b],R[bp]);  
+    std::cout << 0 <<" : R_0 R_0 "<< b <<" "<<bp<<" "<<innerProduct(R[b],R[bp]) <<std::endl;
+  }
+#endif
+
+  GridStopWatch sliceInnerTimer;
+  GridStopWatch sliceMaddTimer;
+  GridStopWatch MatrixTimer;
+  GridStopWatch SolverTimer;
+  SolverTimer.Start();
+
+  int k;
+  int if_print =0;
+  for (k = 1; k <= MaxIterations; k++) {
+
+    RealD rrsum=0;
+    for(int b=0;b<Nblock;b++) rrsum+=real(m_rr(b,b));
+
+    if(PrintInterval && (k%PrintInterval)==0){  
+	if_print=1;
+       std::cout << GridLogMessage << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
+	      <<" / "<<std::sqrt(rrsum/sssum) <<std::endl;
+    } else {
+    if_print=0;
+    std::cout << GridLogIterative << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
+	      <<" / "<<std::sqrt(rrsum/sssum) <<std::endl;
+    }
+
+    MatrixTimer.Start();
+    for(int b=0;b<Nblock;b++) Linop.HermOp(P[b], AP[b]);
+    MatrixTimer.Stop();
+
+    // Alpha
+    sliceInnerTimer.Start();
+ //   sliceInnerProductMatrix(m_pAp,P,AP,Orthog);
+  InnerProductMatrix(m_pAp,P,AP);
+  HermCheck(m_pAp, "P AP",1,if_print);
+  if(if_print) HermCheck(m_pAp, "P AP",0,if_print);
+#if 0
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<Nblock;bp++) {
+    m_pAp(b,bp) = innerProduct(P[b],AP[bp]);  
+    std::cout << k <<" : m_pAp "<< b <<" "<<bp<<" "<<innerProduct(P[b],AP[bp]) <<std::endl;
+  }
+#endif
+    sliceInnerTimer.Stop();
+    m_pAp_inv = m_pAp.inverse();
+  HermCheck(m_pAp_inv, "inv (P AP)",1,if_print);
+  if(if_print) HermCheck(m_pAp_inv, "inv (P AP)",0,if_print);
+if(if_print)
+{
+    m_alpha = m_pAp*m_pAp_inv;
+  for(int b=0;b<Nblock;b++){
+  for(int bp=0;bp<Nblock;bp++) {
+    std::cout << k <<" : pAp*pAp_inv "<< b <<" "<<bp<<" "<<m_alpha(b,bp)<<std::endl;
+  }
+  }
+}
+    m_alpha   = m_pAp_inv * m_rr ; //alpha_k+1 = (P_k+1^t A P_k+1)^-1 (R_k^t R_k)
+
+    // Psi, R update
+    sliceMaddTimer.Start();
+//    sliceMaddMatrix(Psi,m_alpha, P,Psi,Orthog);     // X_k+1=X_k+P_k+1 alpha_k+1
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<Nblock;bp++) {
+    Psi[b] += m_alpha(bp,b)*P[bp];  // X_k+1 = X_k + P_k+1 alpha_k+1
+  }
+
+  for(int b=0;b<Nblock;b++) TMP[b] = R[b];
+//    sliceMaddMatrix(R  ,m_alpha,AP,  R,Orthog,-1.0);// sub alpha * AP to resid
+  for(int b=0;b<Nblock;b++)
+  for(int bp=0;bp<Nblock;bp++) {
+    R[b] -= m_alpha(bp,b)*AP[bp];  // R_k+1 = R_k - AP_k+1 alpha_k+1
+  }
+    sliceMaddTimer.Stop();
+if(if_print)
+{
+//check
+  for(int b=0;b<Nblock;b++){
+  for(int bp=0;bp<Nblock;bp++) {
+    std::cout << k <<" : R_k+1 R_k "<< b <<" "<<bp<<" "<<innerProduct(R[b],TMP[bp]) <<std::endl;
+    std::cout << k <<" : R_k R_k "<< b <<" "<<bp<<" "<<innerProduct(TMP[b],TMP[bp]) <<std::endl;
+  }
+  }
+}
+
+    // Beta
+    m_rr_inv = m_rr.inverse(); //m_rr_inv = (R_k^t R_k)^-1
+    HermCheck(m_rr_inv,"m_rr_inv",1,if_print);
+    if(if_print) HermCheck(m_rr_inv,"m_rr_inv",0,if_print);
+    sliceInnerTimer.Start();
+//    sliceInnerProductMatrix(m_rr,R,R,Orthog);
+    InnerProductMatrix(m_rr,R,R);
+  HermCheck(m_rr,"m_rr",1,if_print);
+  if(if_print) HermCheck(m_rr,"m_rr",0,if_print);
+    sliceInnerTimer.Stop();
+    m_beta = m_rr_inv *m_rr; // beta_k+2 = (R_k^t R_k)^-1 (R_k+1^5 R_k+1)
+//  HermCheck(m_beta,"m_beta");
+
+    // Search update
+    sliceMaddTimer.Start();
+//    sliceMaddMatrix(AP,m_beta,P,R,Orthog);
+  for(int b=0;b<Nblock;b++){
+    AP[b] = R[b];
+  for(int bp=0;bp<Nblock;bp++) {
+    AP[b] += m_beta(bp,b)*P[bp]; //AP = R_k+1 + P_k+1 beta_k+1
+  }
+  }
+if(if_print)
+{
+//check
+    for(int b=0;b<Nblock;b++) Linop.HermOp(P[b], TMP[b]);
+  for(int b=0;b<Nblock;b++){
+  for(int bp=0;bp<Nblock;bp++) {
+    std::cout << k <<" : P_k+2 A P "<< b <<" "<<bp<<" "<<innerProduct(AP[b],TMP[bp]) <<std::endl;
+  }
+  }
+}
+    sliceMaddTimer.Stop();
+  for(int b=0;b<Nblock;b++) P[b]= AP[b]; //P_k+2 = AP
+
+    /*********************
+     * convergence monitor
+     *********************
+     */
+    RealD max_resid=0;
+    RealD rr;
+    for(int b=0;b<Nblock;b++){
+      rr = real(m_rr(b,b))/ssq[b];
+      if ( rr > max_resid ) max_resid = rr;
+    }
+    
+    if ( max_resid < Tolerance*Tolerance ) { 
+
+      SolverTimer.Stop();
+
+      std::cout << GridLogMessage<<"BlockCG converged in "<<k<<" iterations"<<std::endl;
+      for(int b=0;b<Nblock;b++){
+	std::cout << GridLogMessage<< "\t\tblock "<<b<<" computed resid "
+		  << std::sqrt(real(m_rr(b,b))/ssq[b])<<std::endl;
+      }
+	      std::cout << GridLogMessage<<"\tMax residual is "<<std::sqrt(max_resid)<<std::endl;
+
+	  for(int b=0;b<Nblock;b++) { 
+	      Linop.HermOp(Psi[b], AP[b]);
+	AP[b] = AP[b]-Src[b];
+	      std::cout << GridLogMessage <<"\t True residual is " << b<<" "<<std::sqrt(norm2(AP[b])/norm2(Src[b])) <<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 << "\tInnerProd  " << sliceInnerTimer.Elapsed() <<std::endl;
+      std::cout << GridLogMessage << "\tMaddMatrix " << sliceMaddTimer.Elapsed()  <<std::endl;
+	    
+      IterationsToComplete = k;
+      return;
+    }
+
+  }
+  std::cout << GridLogMessage << "BlockConjugateGradient did NOT converge" << std::endl;
+
+  if (ErrorOnNoConverge) assert(0);
+  IterationsToComplete = k;
+}
+
 };
 
 }

lib/algorithms/iterative/InexactPrecConjugateGradient.h

@@ -1,278 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/algorithms/iterative/InexactPrecConjugateGradient.h
-
-    Copyright (C) 2015
-
-Author: Christopher Kelly <ckelly@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_INEXACT_PREC_CONJUGATE_GRADIENT_H_
-#define GRID_INEXACT_PREC_CONJUGATE_GRADIENT_H_
-
-namespace Grid {
-
-//Inexact preconditioned CG based on Golub, Ye, SIAM J. Sci. Comput., 21(4), 1305–1320. 
-//(https://pdfs.semanticscholar.org/d2a9/d5bab02146a7fe3a244677432d21e33a2d98.pdf)
-template <class Field>
-class InexactPreconditionedConjugateGradient : public OperatorFunction<Field> {
- public:
-  bool ErrorOnNoConverge;  // throw an assert when the CG fails to converge.
-                           // Defaults true.
-  RealD Tolerance;
-  Integer MaxIterations;
-  Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion
-  
-  LinearOperatorBase<Field> &Prec;
-
-  InexactPreconditionedConjugateGradient(LinearOperatorBase<Field> &_Prec, RealD tol, Integer maxit, bool err_on_no_conv = true)
-    : Prec(_Prec),
-    Tolerance(tol),
-    MaxIterations(maxit),
-    ErrorOnNoConverge(err_on_no_conv){};
-
-  void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
-
-    psi.checkerboard = src.checkerboard;
-    conformable(psi, src);
-
-    Real ssq = norm2(src);
-    RealD rsq = Tolerance * Tolerance * ssq; //inner stopping condition
-
-    Field p(src);
-    Field r(src);
-    Field rnm1(src);
-    Field mmp(src);
-    Field z(src);
-
-    //Initialize variables
-    Linop.HermOp(psi, mmp);   
-    r = src - mmp;
-    
-    Real cp = norm2(r);
-    
-    p = zero;
-    Real alpha = 0, beta = 0;
-    
-    Real z_nm1_dot_r_nm1;
-
-    int n;
-    for(n=1; n <= MaxIterations; n++) {
-      //Check stopping condition
-      if (cp <= rsq) {
-        Linop.HermOp(psi, mmp);
-        r = mmp - src;
-
-        RealD srcnorm = sqrt(norm2(src));
-        RealD resnorm = sqrt(norm2(r));
-        RealD true_residual = resnorm / srcnorm;
-
-        std::cout << GridLogMessage << "InexactPreconditionedConjugateGradient Converged on iteration " << n << std::endl;
-        std::cout << GridLogMessage << "\tComputed residual " << sqrt(cp / ssq)<<std::endl;
-  	std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
-  	std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
-
-        if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
-
-  	IterationsToComplete = n;	
-
-        return;
-      }
-      
-      std::cout << GridLogIterative << std::setprecision(8)
-		<< "InexactPreconditionedConjugateGradient: n=" << n << " residual " << cp << " target " << rsq << std::endl;
-
-      //Apply preconditioner to current residual
-      Prec.HermOp(r, z);
-
-      //Update beta and store appropriate variables for next iteration
-      Real z_n_dot_r_n  = sqrt(norm(innerProduct(z,r)));
-
-      if(n>1){
-	//  z^T_n ( r_n - r_{n-1} )
-	//  -----------------------
-	//     z^T_{n-1} r_{n-1}
-
-	Real z_n_dot_r_nm1 = sqrt(norm(innerProduct(z,rnm1)));
-	beta = ( z_n_dot_r_n - z_n_dot_r_nm1 ) / z_nm1_dot_r_nm1;
-	std::cout << GridLogIterative << "beta " << beta << std::endl;
-      }
-
-      z_nm1_dot_r_nm1 = z_n_dot_r_n; //for next iteration
-      rnm1 = r;
-
-      axpy(p, beta, p, z); //p = beta * p + z
-
-      //Compute alpha
-      Linop.HermOp(p, mmp);      
-      alpha = z_n_dot_r_n / sqrt(norm(innerProduct(p, mmp)));
-
-      std::cout << GridLogIterative << "alpha " << alpha << std::endl;
-
-      //Update residual and solution
-      cp = axpy_norm(r, -alpha, mmp, r);
-      axpy(psi, alpha, p, psi);
-    }
-    std::cout << GridLogMessage << "InexactPreconditionedConjugateGradient did NOT converge"
-              << std::endl;
-
-    if (ErrorOnNoConverge) assert(0);
-    IterationsToComplete = n;
-  }
-};
-
-
-template<class Field>
-class PolynomialPreconditioner :  public LinearOperatorBase<Field> {
-  Chebyshev<Field> Cheby;
-  LinearOperatorBase<Field> &linop;
-public:
-  int InnerIterations;
-  int order;
-  PolynomialPreconditioner(LinearOperatorBase<Field> &_linop,RealD lo, RealD hi, int _order) 
-    : linop(_linop), Cheby(lo,hi,_order,__InverseApproximation) 
-    {
-      InnerIterations=0;
-      order = _order;
-    };
-
-  void OpDiag (const Field &in, Field &out){ assert(0); }
-  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
-  void Op     (const Field &in, Field &out){ assert(0); }
-  void AdjOp  (const Field &in, Field &out){ assert(0); }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ 
-    HermOp(in,out);
-    n1 = 0; n2 = norm2(out); 
-  }
-  void HermOp(const Field &in, Field &out){ 
-    Cheby(linop,in,out); 
-    InnerIterations+=order;
-  }
-};
-
-template<class Field>
-class DoNothingLinearOperator :  public LinearOperatorBase<Field> {
-public:
-  void OpDiag (const Field &in, Field &out){ assert(0); }
-  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
-  void Op     (const Field &in, Field &out){ assert(0); }
-  void AdjOp  (const Field &in, Field &out){ assert(0); }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ out = in; n1 = 0; n2 = norm2(out); }
-  void HermOp(const Field &in, Field &out){ out = in; }
-};
-
-template<class Field>
-class FixedIterConjugateGradientPreconditioner :  public LinearOperatorBase<Field> {
-public:
-  LinearOperatorBase<Field> &linop;
-  ConjugateGradient<Field> CG;
-
-  FixedIterConjugateGradientPreconditioner (LinearOperatorBase<Field> &_linop, Integer _iter): linop(_linop), CG(1e-20, _iter){
-    CG.ErrorOnNoConverge = false;
-  }
-
-  void OpDiag (const Field &in, Field &out){ assert(0); }
-  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
-  void Op     (const Field &in, Field &out){ assert(0); }
-  void AdjOp  (const Field &in, Field &out){ assert(0); }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ 
-    out = zero;
-    CG(linop,in,out);    
-    n2 = norm2(out);
-  }
-  void HermOp(const Field &in, Field &out){
-    out = zero;
-    CG(linop,in,out);    
-  }
-};
-
-template<class Field>
-class SloppyConjugateGradientPreconditioner :  public LinearOperatorBase<Field> {
-public:
-  LinearOperatorBase<Field> &linop;
-  ConjugateGradient<Field> CG;
-  int InnerIterations;
-
-  SloppyConjugateGradientPreconditioner (LinearOperatorBase<Field> &_linop, Real _resid, Integer max_iter): linop(_linop), CG(_resid, max_iter), InnerIterations(0){
-  }
-
-  void ResetCounters(){ InnerIterations = 0; }
-
-  void OpDiag (const Field &in, Field &out){ assert(0); }
-  void OpDir  (const Field &in, Field &out,int dir,int disp){ assert(0); }
-  void Op     (const Field &in, Field &out){ assert(0); }
-  void AdjOp  (const Field &in, Field &out){ assert(0); }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ 
-    out = zero;
-    CG(linop,in,out);    
-    InnerIterations += CG.IterationsToComplete;
-    n2 = norm2(out);
-  }
-  void HermOp(const Field &in, Field &out){
-    out = zero;
-    CG(linop,in,out);
-    InnerIterations += CG.IterationsToComplete;
-  }
-};
-
-
-template<class FieldH, class FieldL>
-class SloppyConjugateGradientLowerPrecPreconditioner :  public LinearOperatorBase<FieldH> {
-public:
-  LinearOperatorBase<FieldL> &linop;
-  ConjugateGradient<FieldL> CG;
-  GridBase* L_grid; //lower-prec Grid
-  int InnerIterations;
-  FieldL tmp_l1;
-  FieldL tmp_l2;
-
-  SloppyConjugateGradientLowerPrecPreconditioner (LinearOperatorBase<FieldL> &_linop, GridBase* _L_grid, Real _resid, Integer max_iter): linop(_linop), CG(_resid, max_iter), InnerIterations(0), L_grid(_L_grid), tmp_l1(_L_grid), tmp_l2(_L_grid){
-    CG.ErrorOnNoConverge = false;
-  }
-
-  void ResetCounters(){ InnerIterations = 0; }
-
-  void OpDiag (const FieldH &in, FieldH &out){ assert(0); }
-  void OpDir  (const FieldH &in, FieldH &out,int dir,int disp){ assert(0); }
-  void Op     (const FieldH &in, FieldH &out){ assert(0); }
-  void AdjOp  (const FieldH &in, FieldH &out){ assert(0); }
-  void HermOpAndNorm(const FieldH &in, FieldH &out,RealD &n1,RealD &n2){ 
-    precisionChange(tmp_l1, in);
-    tmp_l2 = zero;
-    CG(linop,tmp_l1,tmp_l2);    
-    InnerIterations += CG.IterationsToComplete;
-    precisionChange(out, tmp_l2);
-    n2 = norm2(tmp_l2);
-  }
-  void HermOp(const FieldH &in, FieldH &out){
-    precisionChange(tmp_l1, in);
-    tmp_l2 = zero;
-    CG(linop,tmp_l1,tmp_l2);    
-    InnerIterations += CG.IterationsToComplete;
-    precisionChange(out, tmp_l2);
-  }
-};
-
-}
-
-
-
-#endif

lib/cartesian/Cartesian_red_black.h

@@ -179,6 +179,10 @@ public:
       assert(checker_dim_mask.size() == _ndimension);
       assert(processor_grid.size() == _ndimension);
       assert(simd_layout.size() == _ndimension);
+         std::cout <<"dimensions processor_grid simd_layout checker_dim_mask"<<std::endl;
+      for(int i=0;i<_ndimension;i++){
+         std::cout <<i << " "<<dimensions[i]<<" "<<processor_grid[i]<<" "<< simd_layout[i]<<" "<< checker_dim_mask[i]<<std::endl;
+      }
 
       _fdimensions.resize(_ndimension);
       _gdimensions.resize(_ndimension);

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) {}
   };
 

tests/solver/Test_compression.cc

@@ -1,963 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_dwf_cg_prec.cc
-
-    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 */
-#include<bitset>
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-//Preconditioning:   M psi = chi
-//               =  M P^-1 P psi = chi
-//               =  M P^-1 psi' = chi
-
-//Solve for psi' using M P^-1 as operator, then apply  P^-1 psi' = psi 
-
-//Inexact preconditioned CG requires slight modification because we want to avoid computing P^-1 exactly
-
-
-/////////////////////////////////////////////////////////////
-// Base classes for iterative processes based on operators
-// single input vec, single output vec.
-/////////////////////////////////////////////////////////////
-
-//The compressor
-#if 0
-
-
-//Basic copy of WilsonCompressor for demonstration
-template<class _Hspinor,class _Spinor, class projector>
-class WilsonTestCompressorTemplate
-{
- public:
-  
-  int mu,dag;  
-
-  void Point(int p) { mu=p; };
-
-  WilsonTestCompressorTemplate(int _dag=0){
-    //printf("WilsonTestCompressorTemplate constructor\n");
-    dag = _dag;
-  }
-
-  typedef _Spinor         SiteSpinor;
-  typedef _Hspinor     SiteHalfSpinor;
-  typedef _Hspinor SiteHalfCommSpinor;
-  typedef typename SiteHalfSpinor::vector_type     vComplexIn;
-  constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexIn); //number of complex numbers in SiteHalfSpinor
-
-  inline int CommDatumSize(void) {
-    //printf("WilsonTestCompressorTemplate CommDatumSize\n");
-    return sizeof(SiteHalfCommSpinor);
-  }
-
-  /*****************************************************/
-  /* Compress includes precision change if mpi data is not same */
-  /*****************************************************/
-  inline void Compress(SiteHalfSpinor *buf,Integer o,const SiteSpinor &in) {
-    //printf("WilsonTestCompressorTemplate Compress\n");
-    projector::Proj(buf[o],in,mu,dag);
-  }
-
-  /*****************************************************/
-  /* Exchange includes precision change if mpi data is not same */
-  /*****************************************************/
-  inline void Exchange(SiteHalfSpinor *mp,
-                       SiteHalfSpinor *vp0,
-                       SiteHalfSpinor *vp1,
-		       Integer type,Integer o){
-    //printf("WilsonTestCompressorTemplate Exchange\n");
-    exchange(mp[2*o],mp[2*o+1],vp0[o],vp1[o],type);
-  }
-
-  /*****************************************************/
-  /* Have a decompression step if mpi data is not same */
-  /*****************************************************/
-  inline void Decompress(SiteHalfSpinor *out,
-			 SiteHalfSpinor *in, Integer o) {    
-    //printf("WilsonTestCompressorTemplate Decompress\n");
-    assert(0);
-  }
-
-  /*****************************************************/
-  /* Compress Exchange                                 */
-  /*****************************************************/
-  inline void CompressExchange(SiteHalfSpinor *out0,
-			       SiteHalfSpinor *out1,
-			       const SiteSpinor *in,
-			       Integer j,Integer k, Integer m,Integer type){
-    //printf("WilsonTestCompressorTemplate CompressExchange\n");
-    SiteHalfSpinor temp1, temp2,temp3,temp4;
-    projector::Proj(temp1,in[k],mu,dag);
-    projector::Proj(temp2,in[m],mu,dag);
-    exchange(out0[j],out1[j],temp1,temp2,type);
-  }
-
-  /*****************************************************/
-  /* Pass the info to the stencil */
-  /*****************************************************/
-  inline bool DecompressionStep(void) { return false; }
-
-};
-
-#elif 0
-
-
-//Compressor that unpacks vectorized data to scalar
-template<class _Hspinor,class _Spinor, class projector>
-class WilsonTestCompressorTemplate
-{
- public:
-  
-  int mu,dag;  
-
-  void Point(int p) { mu=p; };
-
-  WilsonTestCompressorTemplate(int _dag=0){
-    dag = _dag;
-  }
-
-  typedef _Spinor         SiteSpinor;
-  typedef _Hspinor     SiteHalfSpinor;
-  typedef _Hspinor SiteHalfCommSpinor;
-  typedef typename SiteHalfSpinor::vector_type     vComplexIn;
-  constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexIn); //number of complex numbers in SiteHalfSpinor
-
-  typedef typename SiteHalfSpinor::scalar_object ScalarSiteHalfSpinor;
-
-  constexpr static int Nsimd = vComplexIn::Nsimd();
-
-  inline int CommDatumSize(void) {
-    return Nsimd*sizeof(ScalarSiteHalfSpinor);
-  }
-
-  /*****************************************************/
-  /* Compress includes precision change if mpi data is not same */
-  /*****************************************************/
-  inline void Compress(SiteHalfSpinor *buf,Integer o,const SiteSpinor &in) {
-    SiteHalfSpinor hsp;
-    projector::Proj(hsp,in,mu,dag);
-    
-    ScalarSiteHalfSpinor* to = (ScalarSiteHalfSpinor*)buf + o*Nsimd;
-    
-    std::vector<ScalarSiteHalfSpinor*> extract_args(Nsimd);
-    for(int i=0;i<Nsimd;i++) extract_args[i] = to+i;
-    extract1(hsp,extract_args,0);
-  }
-
-  /*****************************************************/
-  /* Exchange includes precision change if mpi data is not same */
-  /*****************************************************/
-  inline void Exchange(SiteHalfSpinor *mp,
-                       SiteHalfSpinor *vp0,
-                       SiteHalfSpinor *vp1,
-		       Integer type,Integer o){
-    ScalarSiteHalfSpinor* vpp0 = (ScalarSiteHalfSpinor*)vp0 + o*Nsimd;
-    ScalarSiteHalfSpinor* vpp1 = (ScalarSiteHalfSpinor*)vp1 + o*Nsimd;
-    
-    std::vector<ScalarSiteHalfSpinor*> merge_args0(Nsimd), merge_args1(Nsimd);
-    for(int i=0;i<Nsimd;i++){
-      merge_args0[i] = vpp0+i;
-      merge_args1[i] = vpp1+i;
-    }
-
-    SiteHalfSpinor vt0,vt1;
-    merge1(vt0,merge_args0,0);
-    merge1(vt1,merge_args1,0);
-
-    exchange(mp[2*o],mp[2*o+1],vt0,vt1,type);
-  }
-
-  /*****************************************************/
-  /* Have a decompression step if mpi data is not same */
-  /*****************************************************/
-  inline void Decompress(SiteHalfSpinor *out,
-			 SiteHalfSpinor *in, Integer o) {    
-    ScalarSiteHalfSpinor* hin = (ScalarSiteHalfSpinor*)in + o*Nsimd;
-    std::vector<ScalarSiteHalfSpinor*> merge_args(Nsimd);
-    for(int i=0;i<Nsimd;i++) merge_args[i] = hin+i;
-    merge1(out[o],merge_args,0);
-  }
-
-  /*****************************************************/
-  /* Compress Exchange                                 */
-  /*****************************************************/
-  inline void CompressExchange(SiteHalfSpinor *out0,
-			       SiteHalfSpinor *out1,
-			       const SiteSpinor *in,
-			       Integer j,Integer k, Integer m,Integer type){
-    SiteHalfSpinor temp1, temp2,temp3,temp4;
-    projector::Proj(temp1,in[k],mu,dag);
-    projector::Proj(temp2,in[m],mu,dag);
-    exchange(temp3,temp4,temp1,temp2,type);
-
-    ScalarSiteHalfSpinor* hout0 = (ScalarSiteHalfSpinor*)out0 + j*Nsimd;
-    ScalarSiteHalfSpinor* hout1 = (ScalarSiteHalfSpinor*)out1 + j*Nsimd;
-
-    std::vector<ScalarSiteHalfSpinor*> extract_args0(Nsimd), extract_args1(Nsimd);
-    for(int i=0;i<Nsimd;i++){
-      extract_args0[i] = hout0+i;
-      extract_args1[i] = hout1+i;
-    }
-    extract1(temp3,extract_args0,0);
-    extract1(temp4,extract_args1,0);
-  }
-
-  /*****************************************************/
-  /* Pass the info to the stencil */
-  /*****************************************************/
-  inline bool DecompressionStep(void) { return true; }
-
-};
-
-#else
-
-//Access elements of std::complex
-template<typename T>
-inline T & cmplx_reim(std::complex<T> &c, const int reim){
-  return reinterpret_cast<T(&)[2]>(c)[reim];
-}
-
-template<typename T>
-inline const T & cmplx_reim(const std::complex<T> &c, const int reim){
-  return reinterpret_cast<const T(&)[2]>(c)[reim];
-}
-
-
-//Pack and unpack float/double to fixed point representation of SZ bits
-template<int SZ>
-struct signedIntMap{};
-
-template<>
-struct signedIntMap<8>{ typedef int8_t type; };
-template<>
-struct signedIntMap<16>{ typedef int16_t type; };
-
-
-template<typename T, int SZ>
-inline typename signedIntMap<SZ>::type packN(T val){
-  return typename signedIntMap<SZ>::type( (1<<(SZ-2) ) * val );
-}
-template<typename T, int SZ>
-inline T unpackN(typename signedIntMap<SZ>::type val){
-  return T(val)/(1<<(SZ-2));
-}
-
-template<typename T>
-struct getHalfSpinorColors{
-  //template <typename vtype> using iImplHalfSpinor        = iScalar<iVector<iVector<vtype, Dimension>, Nhs> >;
-  enum { value = sizeof(typename T::element::element)/sizeof(typename T::element::element::element) };
-};
-
-//Compressor that compresses to a single magnitude and Nhs*Dimension fixed point integers of size packSize bits
-template<class _Hspinor,class _Spinor, class projector, int packSize = 16>
-class WilsonTestCompressorTemplate
-{
- public:
-  
-  int mu,dag;  
-
-  void Point(int p) { mu=p; };
-
-  WilsonTestCompressorTemplate(int _dag=0){
-    dag = _dag;
-  }
-
-  typedef _Spinor         SiteSpinor;
-  typedef _Hspinor     SiteHalfSpinor;
-  typedef _Hspinor SiteHalfCommSpinor;
-  typedef typename SiteHalfSpinor::vector_type     vComplexIn;
-  constexpr static int Nw=sizeof(SiteHalfSpinor)/sizeof(vComplexIn); //number of complex numbers in SiteHalfSpinor
-
-  typedef typename SiteHalfSpinor::scalar_object ScalarSiteHalfSpinor;
-
-  constexpr static int Nsimd = vComplexIn::Nsimd();
-  constexpr static int Dimension = getHalfSpinorColors<SiteHalfSpinor>::value;
- 
-  typedef typename ScalarSiteHalfSpinor::scalar_type stype; //std::complex
-  typedef typename stype::value_type srtype; //float/double
-
-  //Pack and unpack *scalar* SiteHalfSpinor objects
-  void packSpinor(void* tov, const ScalarSiteHalfSpinor &from){
-    uint8_t* to = (uint8_t*)tov;
-    typedef typename signedIntMap<packSize>::type packedType;
-
-    srtype max = 0;
-    for(int s=0;s<Nhs;s++)
-      for(int c=0;c<Dimension;c++)
-	for(int reim=0;reim<2;reim++)
-	  if(fabs(cmplx_reim( from()(s)(c), reim )) > max )
-	    max =  fabs(cmplx_reim( from()(s)(c), reim )) ;
-  
-    *( (srtype*)to ) = max; //copy the normalization to the buffer
-    to += sizeof(srtype);
-  
-    packedType *top = (packedType*)to;
-    packedType p;
-    srtype q;
-    for(int s=0;s<Nhs;s++)
-      for(int c=0;c<Dimension;c++)
-	for(int reim=0;reim<2;reim++){
-	  q = cmplx_reim( from()(s)(c), reim );
-	  if(max != 0.) q /= max;
-	  *(top++) = packN<srtype,packSize>(q);
-	}
-  }
-
-  void packSpinor(void* tov, const SiteHalfSpinor &from){
-    uint8_t* to = (uint8_t*)tov;
-    std::vector<ScalarSiteHalfSpinor> extracted(Nsimd);
-    extract(from,extracted);
-
-    static const int incr = sizeof(srtype) + Nhs*Dimension*2*sizeof(typename signedIntMap<packSize>::type);
-
-    for(int i=0;i<Nsimd;i++){
-      packSpinor((void*)to, extracted[i]);
-      to += incr;
-    }
-  }
-
-
-  void unpackSpinor(ScalarSiteHalfSpinor &to, void* fromv){
-    uint8_t* from = (uint8_t*)fromv;
-    typedef typename signedIntMap<packSize>::type packedType;
-
-    srtype norm = *( (srtype*)from ); 
-    from += sizeof(srtype);
-
-    packedType *fromp = (packedType*)from;
-    srtype q;
-    for(int s=0;s<Nhs;s++)
-      for(int c=0;c<Dimension;c++)
-	for(int reim=0;reim<2;reim++){
-	  q = unpackN<srtype,packSize>(*(fromp++) );
-	  if(norm != 0.) q *= norm;
-	  cmplx_reim( to()(s)(c), reim ) = q;
-	}
-  }
-
-  void unpackSpinor(SiteHalfSpinor &to, void* fromv){
-    uint8_t* from = (uint8_t*)fromv;
-    std::vector<ScalarSiteHalfSpinor> unpacked(Nsimd);
-
-    static const int incr = sizeof(srtype) + Nhs*Dimension*2*sizeof(typename signedIntMap<packSize>::type);
-
-    for(int i=0;i<Nsimd;i++){
-      unpackSpinor(unpacked[i],(void*)from);
-      from += incr;
-    }
-
-    merge(to,unpacked);
-  }
-
-  inline int CommDatumSize(void) {
-    return Nsimd*(  sizeof(srtype) + Nhs*Dimension*2*sizeof(typename signedIntMap<packSize>::type) );
-  }
-
-  /*****************************************************/
-  /* Compress includes precision change if mpi data is not same */
-  /*****************************************************/
-  void Compress(SiteHalfSpinor *buf,Integer o,const SiteSpinor &in) {
-    SiteHalfSpinor hsp;
-    projector::Proj(hsp,in,mu,dag);
-    
-    uint8_t* to = (uint8_t*)buf + o*CommDatumSize();
-    packSpinor(to, hsp);
-  }
-
-  /*****************************************************/
-  /* Exchange includes precision change if mpi data is not same */
-  /*****************************************************/
-  void Exchange(SiteHalfSpinor *mp,
-                       SiteHalfSpinor *vp0,
-                       SiteHalfSpinor *vp1,
-		       Integer type,Integer o){
-    uint8_t* vpp0 = (uint8_t*)vp0 + o*CommDatumSize();
-    uint8_t* vpp1 = (uint8_t*)vp1 + o*CommDatumSize();
-
-    SiteHalfSpinor vt0, vt1;
-    unpackSpinor(vt0, vpp0);
-    unpackSpinor(vt1, vpp1);
-
-    exchange(mp[2*o],mp[2*o+1],vt0,vt1,type);
-  }
-
-  /*****************************************************/
-  /* Have a decompression step if mpi data is not same */
-  /*****************************************************/
-  void Decompress(SiteHalfSpinor *out,
-			 SiteHalfSpinor *in, Integer o) {    
-    uint8_t* hin = (uint8_t*)in + o*CommDatumSize();
-    unpackSpinor(out[o],hin);
-  }
-
-  /*****************************************************/
-  /* Compress Exchange                                 */
-  /*****************************************************/
-  void CompressExchange(SiteHalfSpinor *out0,
-			       SiteHalfSpinor *out1,
-			       const SiteSpinor *in,
-			       Integer j,Integer k, Integer m,Integer type){
-    SiteHalfSpinor temp1, temp2,temp3,temp4;
-    projector::Proj(temp1,in[k],mu,dag);
-    projector::Proj(temp2,in[m],mu,dag);
-    exchange(temp3,temp4,temp1,temp2,type);
-
-    uint8_t* hout0 = (uint8_t*)out0 + j*CommDatumSize();
-    uint8_t* hout1 = (uint8_t*)out1 + j*CommDatumSize();
-    packSpinor(hout0, temp3);
-    packSpinor(hout1, temp4);
-  }
-
-  /*****************************************************/
-  /* Pass the info to the stencil */
-  /*****************************************************/
-  inline bool DecompressionStep(void) { return true; }
-
-};
-
-
-#endif
-
-
-
-template<typename HS,typename S, int packSize> using WilsonTestCompressor = WilsonTestCompressorTemplate<HS,S,WilsonProjector,packSize>;
-
-
-template<class vobj,class cobj>
-class WilsonStencilBasic : public CartesianStencil<vobj,cobj> {
-public:
-  double timer0;
-  double timer1;
-  double timer2;
-  double timer3;
-  double timer4;
-  double timer5;
-  double timer6;
-  uint64_t callsi;
-  void ZeroCountersi(void)
-  {
-    timer0=0;
-    timer1=0;
-    timer2=0;
-    timer3=0;
-    timer4=0;
-    timer5=0;
-    timer6=0;
-    callsi=0;
-  }
-  void Reporti(int calls)
-  {
-    if ( timer0 ) std::cout << GridLogMessage << " timer0 (HaloGatherOpt) " <<timer0/calls <<std::endl;
-    if ( timer1 ) std::cout << GridLogMessage << " timer1 (Communicate)   " <<timer1/calls <<std::endl;
-    if ( timer2 ) std::cout << GridLogMessage << " timer2 (CommsMerge )   " <<timer2/calls <<std::endl;
-    if ( timer3 ) std::cout << GridLogMessage << " timer3 (commsMergeShm) " <<timer3/calls <<std::endl;
-    if ( timer4 ) std::cout << GridLogMessage << " timer4 " <<timer4 <<std::endl;
-  }
-
-
-  std::vector<int> same_node;
-  std::vector<int> surface_list;
-
-  WilsonStencilBasic(GridBase *grid,
-		int npoints,
-		int checkerboard,
-		const std::vector<int> &directions,
-		const std::vector<int> &distances)  
-    : CartesianStencil<vobj,cobj> (grid,npoints,checkerboard,directions,distances) ,
-    same_node(npoints)
-  { 
-    ZeroCountersi();
-    surface_list.resize(0);
-  };
-
-  void BuildSurfaceList(int Ls,int vol4){
-
-    // find same node for SHM
-    // Here we know the distance is 1 for WilsonStencil
-    for(int point=0;point<this->_npoints;point++){
-      same_node[point] = this->SameNode(point);
-    }
-    
-    for(int site = 0 ;site< vol4;site++){
-      int local = 1;
-      for(int point=0;point<this->_npoints;point++){
-	if( (!this->GetNodeLocal(site*Ls,point)) && (!same_node[point]) ){ 
-	  local = 0;
-	}
-      }
-      if(local == 0) { 
-	surface_list.push_back(site);
-      }
-    }
-  }
-
-
-  template < class compressor>
-  void HaloExchangeOpt(const Lattice<vobj> &source,compressor &compress) 
-  {
-    std::vector<std::vector<CommsRequest_t> > reqs;
-    this->HaloExchangeOptGather(source,compress);
-    double t1=usecond();
-    this->Communicate();
-    double t2=usecond(); timer1 += t2-t1;
-    this->CommsMerge(compress);
-    double t3=usecond(); timer2 += t3-t2;
-    this->CommsMergeSHM(compress);
-    double t4=usecond(); timer3 += t4-t3;
-  }
-  
-  template <class compressor>
-  void HaloExchangeOptGather(const Lattice<vobj> &source,compressor &compress){
-    this->Prepare();
-    double t0=usecond();
-    this->HaloGatherOpt(source,compress);
-    double t1=usecond();
-    timer0 += t1-t0;
-    callsi++;
-  }
-
-  template <class compressor>
-  void HaloGatherOpt(const Lattice<vobj> &source,compressor &compress)
-  {
-    this->halogtime-=usecond();
-    this->HaloGather(source,compress);
-    this->halogtime+=usecond();
-  }
-};
-
-
-
-
-//This is hideous
-template<class S, int packSize = 16>
-class WilsonCompressedCommsImpl: public WilsonImpl<S,FundamentalRepresentation,CoeffReal>{
-public:
-  typedef WilsonImpl<S,FundamentalRepresentation,CoeffReal> WilsonBase;
-
-#define INHERIT_BASE(TYPE) typedef typename WilsonBase::TYPE TYPE
-
-  INHERIT_BASE(Gimpl);
-  INHERIT_GIMPL_TYPES(Gimpl);
-
-  INHERIT_BASE(Coeff_t);
-  
-  INHERIT_BASE(SiteSpinor);
-  INHERIT_BASE(SitePropagator);
-  INHERIT_BASE(SiteHalfSpinor);
-  INHERIT_BASE(SiteHalfCommSpinor);    
-  INHERIT_BASE(SiteDoubledGaugeField);
-
-  INHERIT_BASE(FermionField);
-  INHERIT_BASE(PropagatorField);
-  INHERIT_BASE(DoubledGaugeField);
-
-  //typedef WilsonCompressor<SiteHalfCommSpinor,SiteHalfSpinor, SiteSpinor> Compressor;
-  typedef WilsonTestCompressor<SiteHalfSpinor, SiteSpinor, packSize> Compressor;  
-
-  INHERIT_BASE(ImplParams);
-  //INHERIT_BASE(StencilImpl);
-  typedef WilsonStencilBasic<SiteSpinor, SiteHalfSpinor> StencilImpl;
-
-  WilsonCompressedCommsImpl(const ImplParams &p = ImplParams()) : WilsonBase(p){}
-  
-  inline void multLink(SiteHalfSpinor &phi,
-		       const SiteDoubledGaugeField &U,
-		       const SiteHalfSpinor &chi,
-		       int mu,
-		       StencilEntry *SE,
-		       StencilImpl &St) {
-    mult(&phi(), &U(mu), &chi());
-  }
-
-#undef INHERIT_BASE
-};
-
-typedef WilsonCompressedCommsImpl<vComplexF,8> WilsonCompressedComms8ImplF;
-typedef WilsonCompressedCommsImpl<vComplexD,8> WilsonCompressedComms8ImplD;
-typedef WilsonCompressedCommsImpl<vComplexF,16> WilsonCompressedComms16ImplF;
-typedef WilsonCompressedCommsImpl<vComplexD,16> WilsonCompressedComms16ImplD;
-
-
-#define TO_INSTANTIATE \
-  DOIT(WilsonCompressedComms8ImplF)\
-  DOIT(WilsonCompressedComms8ImplD)\
-  DOIT(WilsonCompressedComms16ImplF)\
-  DOIT(WilsonCompressedComms16ImplD)
-
-#include "InstantiateImpl.impl"
-
-#undef TO_INSTANTIATE
-
-typedef DomainWallFermion<WilsonCompressedComms8ImplD> DomainWallFermionCompressedComms8D;
-typedef DomainWallFermion<WilsonCompressedComms8ImplF> DomainWallFermionCompressedComms8F;
-typedef DomainWallFermion<WilsonCompressedComms16ImplD> DomainWallFermionCompressedComms16D;
-typedef DomainWallFermion<WilsonCompressedComms16ImplF> DomainWallFermionCompressedComms16F;
-
-template<typename T>
-T parse(const std::string &name, std::istream &in){
-  std::string p;
-  in >> p;
-  assert(p==name);
-  char eq;
-  in >> eq;
-  assert(eq == '=');
-  T out;
-  in >> out;
-  return out;
-}
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  int Ls=8;
-  RealD mass=0.1;
-  RealD outer_tol = 1e-8;
-  RealD inner_tol_full = 1e-5;
-  RealD inner_tol_half = 1e-5;
-  RealD inner_tol_16c = 1e-5;
-  RealD inner_tol_8c = 1e-5;
-
-  RealD relup_delta_full = 0.1;
-  RealD relup_delta_half = 0.1;
-  RealD relup_delta_16c = 0.1;
-  RealD relup_delta_8c = 0.1;
-
-  std::string config_file = "";
-
-  for(int i=1;i<argc;i++){
-    if(std::string(argv[i]) == "--params"){
-      std::ifstream f(argv[i+1]);
-      f.exceptions ( std::ifstream::failbit | std::ifstream::badbit );
-      Ls = parse<int>("Ls", f);
-#define PARSEIT(NM) NM = parse<RealD>(#NM, f)
-      PARSEIT(mass);
-      PARSEIT(outer_tol);
-      PARSEIT(inner_tol_full);
-      PARSEIT(inner_tol_half);
-      PARSEIT(inner_tol_16c);
-      PARSEIT(inner_tol_8c);
-      PARSEIT(relup_delta_full);
-      PARSEIT(relup_delta_half);
-      PARSEIT(relup_delta_16c);
-      PARSEIT(relup_delta_8c);
-#undef PARSEIT
-
-      //f >> outer_tol >> inner_tol_full >> inner_tol_half >> inner_tol_16c >> inner_tol_8c;      
-    }else if(std::string(argv[i]) == "--config"){
-      config_file = argv[i+1];
-    }
-  }
-  
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  GridCartesian         * UGrid_f   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
-  GridCartesian         * FGrid_f   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f);
-  GridRedBlackCartesian * FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f);
-  
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
-
-  LatticeFermionD    src(FGrid); random(RNG5,src);
-  LatticeFermionD result(FGrid); result=zero;
-  LatticeGaugeFieldD Umu(UGrid);
-  LatticeGaugeFieldF Umu_f(UGrid_f); 
-
-  if(config_file.size() > 0){
-    FieldMetaData header;
-    NerscIO::readConfiguration(Umu,header,config_file);
-  }else{
-    SU3::HotConfiguration(RNG4,Umu);
-  }
-
-  precisionChange(Umu_f,Umu);
-  
-  RealD M5=1.8;
-
-  LatticeFermionD    src_o(FrbGrid);
-  pickCheckerboard(Odd,src_o,src);
-
-  if(0){ //Test preconditioned CG
-    LatticeFermionD result_o(FrbGrid);
-    LatticeFermionD result_o_2(FrbGrid);
-    result_o.checkerboard = Odd;
-    result_o = zero;
-    result_o_2.checkerboard = Odd;
-    result_o_2 = zero;
-
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-    //DoNothingLinearOperator<LatticeFermionD> Prec;
-    //FixedIterConjugateGradientPreconditioner<LatticeFermionD> Prec(HermOpEO, 20);
-    SloppyConjugateGradientPreconditioner<LatticeFermionD> Prec(HermOpEO, 1e-2, 1000);
-
-    std::cout << "Preconditioned CG" << std::endl;
-    InexactPreconditionedConjugateGradient<LatticeFermionD> pCG(Prec,1.0e-8,10000);
-    pCG(HermOpEO,src_o,result_o);
-
-    std::cout << "Starting regular CG" << std::endl;
-    ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-    CG(HermOpEO,src_o,result_o_2);
-
-    LatticeFermionD diff_o(FrbGrid);
-    RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
-
-    std::cout << "pCG HermOp applications " << pCG.IterationsToComplete << "(outer) + " << Prec.InnerIterations << "(inner) = " << pCG.IterationsToComplete + Prec.InnerIterations << std::endl;
-    std::cout << "CG HermOp applications " << CG.IterationsToComplete << std::endl;
-    std::cout << "Diff between results: " << diff << std::endl;
-  }
-
-  if(0){ //Test compressor
-    LatticeFermionD result_o(FrbGrid);
-    LatticeFermionD result_o_2(FrbGrid);
-    result_o.checkerboard = Odd;
-    result_o = zero;
-    result_o_2.checkerboard = Odd;
-    result_o_2 = zero;
-
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-
-    DomainWallFermionCompressedComms16D DdwfC(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionCompressedComms16D,LatticeFermionD> HermOpEOC(DdwfC);
-
-    std::cout << "Starting regular CG with compressed operator" << std::endl;
-    Integer iter1;
-    {
-      ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-      CG.ErrorOnNoConverge = false;
-      CG(HermOpEOC,src_o,result_o);
-      iter1 = CG.IterationsToComplete;
-    }
-    Integer iter2;
-    {
-      std::cout << "Starting regular CG" << std::endl;
-      ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-      CG(HermOpEO,src_o,result_o_2);
-      iter2 = CG.IterationsToComplete;
-    }
-
-    LatticeFermionD diff_o(FrbGrid);
-    RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
-
-    std::cout << "CG HermOp CC applications " << iter1 << std::endl;
-    std::cout << "CG HermOp applications " << iter2 << std::endl;
-    std::cout << "Diff between results: " << diff << std::endl;
-  }
-  
-  if(1){ //Compare mixed prec restarted single/single internal with same but with single/compressed
-    LatticeFermionD result_o_full(FrbGrid);
-    LatticeFermionD result_o_half(FrbGrid);
-    LatticeFermionD result_o_16(FrbGrid);
-    LatticeFermionD result_o_8(FrbGrid);
-    result_o_full.checkerboard = Odd;
-    result_o_full = zero;
-    result_o_16 = result_o_8 = result_o_half = result_o_full;
-
-    //Std
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-
-    DomainWallFermionF Ddwf_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionF,LatticeFermionF> HermOpEO_f(Ddwf_f);
-
-    //1/2 prec
-    DomainWallFermionFH Ddwfhalf_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionFH,LatticeFermionF> HermOpEOhalf_f(Ddwfhalf_f);
-
-    //16
-    DomainWallFermionCompressedComms16F DdwfC16_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionCompressedComms16F,LatticeFermionF> HermOpEOC16_f(DdwfC16_f);
-
-    //8
-    DomainWallFermionCompressedComms8F DdwfC8_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionCompressedComms8F,LatticeFermionF> HermOpEOC8_f(DdwfC8_f);
-
-    #define ALGORITHM_MIXEDCG
-    //#define ALGORITHM_RELUP
-    //#define ALGORITHM_SLOPPY_PREC_CG
-
-#ifdef ALGORITHM_MIXEDCG
-    std::cout << "Starting mixed CG with single/compressed-16 inner\n";
-
-    Integer inner_16, outer_16, patchup_16;
-    {
-      MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(outer_tol, 10000, 50, FrbGrid_f, HermOpEOC16_f, HermOpEO);
-      mCG.InnerTolerance = inner_tol_16c;
-      mCG(src_o,result_o_16);
-      inner_16 = mCG.TotalInnerIterations; outer_16 = mCG.TotalOuterIterations; patchup_16 = mCG.TotalFinalStepIterations;
-    }
-
-    std::cout << "Starting mixed CG with single/compressed-8 inner\n";
-    Integer inner_8, outer_8, patchup_8;
-    {
-      MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(outer_tol, 10000, 50, FrbGrid_f, HermOpEOC8_f, HermOpEO);
-      mCG.InnerTolerance = inner_tol_8c;
-      mCG(src_o,result_o_8);
-      inner_8 = mCG.TotalInnerIterations; outer_8 = mCG.TotalOuterIterations; patchup_8 = mCG.TotalFinalStepIterations;
-    }
-
-    std::cout << "Starting mixed CG with single/half inner\n";
-    Integer inner_half, outer_half, patchup_half;
-    {
-      MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(outer_tol, 10000, 50, FrbGrid_f, HermOpEOhalf_f, HermOpEO);
-      mCG.InnerTolerance = inner_tol_half;
-      mCG(src_o,result_o_half);
-      inner_half = mCG.TotalInnerIterations; outer_half = mCG.TotalOuterIterations; patchup_half = mCG.TotalFinalStepIterations;
-    }
-
-    std::cout << "Starting mixed CG with single/single inner\n";
-    Integer inner_full, outer_full, patchup_full;
-    {
-      MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(outer_tol, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO);
-      mCG.InnerTolerance = inner_tol_full;
-      mCG(src_o,result_o_full);
-      inner_full = mCG.TotalInnerIterations; outer_full = mCG.TotalOuterIterations; patchup_full = mCG.TotalFinalStepIterations;
-    }
-#elif defined(ALGORITHM_RELUP)
-    std::cout << "Starting relup CG with single/compressed-16 inner\n";    
-    Integer inner_16, outer_16, patchup_16;
-    {
-      ConjugateGradientReliableUpdate<LatticeFermionD,LatticeFermionF> relup(outer_tol, 2000, relup_delta_16c, FrbGrid_f, HermOpEOC16_f, HermOpEO);
-      relup(src_o,result_o_16);
-      inner_16 = relup.IterationsToComplete; outer_16 = relup.ReliableUpdatesPerformed; patchup_16 = relup.IterationsToCleanup;
-    }
-
-    std::cout << "Starting relup CG with single/compressed-8 inner\n";
-    Integer inner_8, outer_8, patchup_8;
-    {
-      ConjugateGradientReliableUpdate<LatticeFermionD,LatticeFermionF> relup(outer_tol, 2000, relup_delta_8c, FrbGrid_f, HermOpEOC8_f, HermOpEO);
-      relup.ErrorOnNoConverge = false;
-      relup(src_o,result_o_8);
-      inner_8 = relup.IterationsToComplete; outer_8 = relup.ReliableUpdatesPerformed; patchup_8 = relup.IterationsToCleanup;     
-    }
-
-    std::cout << "Starting relup CG with single/half inner\n";
-    Integer inner_half, outer_half, patchup_half;
-    {
-      ConjugateGradientReliableUpdate<LatticeFermionD,LatticeFermionF> relup(outer_tol, 2000, relup_delta_half, FrbGrid_f, HermOpEOhalf_f, HermOpEO);
-      relup(src_o,result_o_half);
-      inner_half = relup.IterationsToComplete; outer_half = relup.ReliableUpdatesPerformed; patchup_half = relup.IterationsToCleanup;
-    }
-
-    std::cout << "Starting relup CG with single/single inner\n";
-    Integer inner_full, outer_full, patchup_full;
-    {
-      ConjugateGradientReliableUpdate<LatticeFermionD,LatticeFermionF> relup(outer_tol, 2000, relup_delta_full, FrbGrid_f, HermOpEO_f, HermOpEO);
-      relup(src_o,result_o_full);
-      inner_full = relup.IterationsToComplete; outer_full = relup.ReliableUpdatesPerformed; patchup_full = relup.IterationsToCleanup;
-    }
-#elif defined(ALGORITHM_SLOPPY_PREC_CG)
-
-    std::cout << "Starting sloppy pCG with single/compressed-16 inner\n";    
-    Integer inner_16, outer_16;
-    {
-      SloppyConjugateGradientLowerPrecPreconditioner<LatticeFermionD,LatticeFermionF> prec(HermOpEOC16_f, FrbGrid_f, inner_tol_16c, 1000);
-      InexactPreconditionedConjugateGradient<LatticeFermionD> CG(prec, outer_tol, 100);
-      CG(HermOpEO,src_o,result_o_16);
-      inner_16 = prec.InnerIterations; outer_16 = CG.IterationsToComplete;
-    }
-
-    std::cout << "Starting sloppy pCG with single/compressed-8 inner\n";    
-    Integer inner_8, outer_8;
-    {
-      SloppyConjugateGradientLowerPrecPreconditioner<LatticeFermionD,LatticeFermionF> prec(HermOpEOC8_f, FrbGrid_f, inner_tol_8c, 1000);
-      InexactPreconditionedConjugateGradient<LatticeFermionD> CG(prec, outer_tol, 100);
-      CG(HermOpEO,src_o,result_o_8);
-      inner_8 = prec.InnerIterations; outer_8 = CG.IterationsToComplete;
-    }
-
-    std::cout << "Starting sloppy pCG with single/half inner\n";    
-    Integer inner_half, outer_half;
-    {
-      SloppyConjugateGradientLowerPrecPreconditioner<LatticeFermionD,LatticeFermionF> prec(HermOpEOhalf_f, FrbGrid_f, inner_tol_half, 1000);
-      InexactPreconditionedConjugateGradient<LatticeFermionD> CG(prec, outer_tol, 100);
-      CG(HermOpEO,src_o,result_o_half);
-      inner_half = prec.InnerIterations; outer_half = CG.IterationsToComplete;
-    }
-
-    std::cout << "Starting sloppy pCG with single/single inner\n";    
-    Integer inner_full, outer_full;
-    {
-      SloppyConjugateGradientLowerPrecPreconditioner<LatticeFermionD,LatticeFermionF> prec(HermOpEO_f, FrbGrid_f, inner_tol_full, 1000);
-      InexactPreconditionedConjugateGradient<LatticeFermionD> CG(prec, outer_tol, 100);
-      CG(HermOpEO,src_o,result_o_full);
-      inner_full = prec.InnerIterations; outer_full = CG.IterationsToComplete;
-    }
-
-#endif
-
-    std::cout << "Ls " << Ls << std::endl;
-    std::cout << "Mass " << mass << std::endl;
-    std::cout << "Outer tolerance " << outer_tol << std::endl;
-
-#if defined(ALGORITHM_MIXEDCG) || defined(ALGORITHM_SLOPPY_PREC_CG)
-    std::cout << "Inner tol full " << inner_tol_full << std::endl;
-    std::cout << "Inner tol 1/2 prec " << inner_tol_half << std::endl;
-    std::cout << "Inner tol compressed-16 " << inner_tol_16c << std::endl;
-    std::cout << "Inner tol compressed-8 " << inner_tol_8c << std::endl;
-#elif defined(ALGORITHM_RELUP)
-    std::cout << "Relup delta full " << relup_delta_full << std::endl;
-    std::cout << "Relup delta 1/2 prec " << relup_delta_half << std::endl;
-    std::cout << "Relup delta compressed-16 " << relup_delta_16c << std::endl;
-    std::cout << "Relup delta compressed-8 " << relup_delta_8c << std::endl;
-#endif
-    LatticeFermionD diff_o(FrbGrid);
-    RealD diff = axpy_norm(diff_o, -1.0, result_o_16, result_o_full);
-    std::cout << "Diff between results (s/c16): " << diff << std::endl;
-
-    diff = axpy_norm(diff_o, -1.0, result_o_8, result_o_full);
-    std::cout << "Diff between results (s/c8): " << diff << std::endl;
-
-    diff = axpy_norm(diff_o, -1.0, result_o_half, result_o_full);
-    std::cout << "Diff between results (s/h): " << diff << std::endl;
-
-#if defined(ALGORITHM_MIXEDCG) || defined(ALGORITHM_RELUP)
-    std::cout << "Iterations (s/c16) inner: " << inner_16 << " outer: " << outer_16 << " patchup: " << patchup_16 << std::endl;
-    std::cout << "Iterations (s/c8) inner: " << inner_8 << " outer: " << outer_8 << " patchup: " << patchup_8 << std::endl;
-    std::cout << "Iterations (s/h) inner: " << inner_half << " outer: " << outer_half << " patchup: " << patchup_half << std::endl;
-    std::cout << "Iterations (s/s) inner: " << inner_full << " outer: " << outer_full << " patchup: " << patchup_full << std::endl;
-#else
-    std::cout << "Iterations (s/c16) inner: " << inner_16 << " outer: " << outer_16 << std::endl;
-    std::cout << "Iterations (s/c8) inner: " << inner_8 << " outer: " << outer_8 << std::endl;
-    std::cout << "Iterations (s/h) inner: " << inner_half << " outer: " << outer_half << std::endl;
-    std::cout << "Iterations (s/s) inner: " << inner_full << " outer: " << outer_full << std::endl;
-#endif
-
-  }
-
-
-  Grid_finalize();
-}

tests/solver/Test_dwf_mrhs_cg.cc

@@ -68,6 +68,7 @@ int main (int argc, char ** argv)
 
   int nrhs = 1;
   int me;
+  for(int i=0;i<mpi_layout.size();i++) cout <<" node split = "<<mpi_layout[i]<<" "<<mpi_split[i]<<endl;
   for(int i=0;i<mpi_layout.size();i++) nrhs *= (mpi_layout[i]/mpi_split[i]);
 
   GridCartesian         * SGrid = new GridCartesian(GridDefaultLatt(),
@@ -99,12 +100,6 @@ int main (int argc, char ** argv)
   // Bounce these fields to disk
   ///////////////////////////////////////////////////////////////
 
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::cout << GridLogMessage << " Writing out in parallel view "<<std::endl;
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  emptyUserRecord record;
-  std::string file("./scratch.scidac");
-  std::string filef("./scratch.scidac.ferm");
 
   LatticeGaugeField s_Umu(SGrid);
   FermionField s_src(SFGrid);
@@ -114,57 +109,10 @@ int main (int argc, char ** argv)
 
   {
     FGrid->Barrier();
-    ScidacWriter _ScidacWriter(FGrid->IsBoss());
-    _ScidacWriter.open(file);
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    std::cout << GridLogMessage << " Writing out gauge field "<<std::endl;
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    _ScidacWriter.writeScidacFieldRecord(Umu,record);
-    _ScidacWriter.close();
-    FGrid->Barrier();
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    std::cout << GridLogMessage << " Reading in gauge field "<<std::endl;
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    ScidacReader  _ScidacReader;
-    _ScidacReader.open(file);
-    _ScidacReader.readScidacFieldRecord(s_Umu,record);
-    _ScidacReader.close();
-    FGrid->Barrier();
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    std::cout << GridLogMessage << " Read in gauge field "<<std::endl;
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
   }
 
 
   {
-    for(int n=0;n<nrhs;n++){
-
-      std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-      std::cout << GridLogMessage << " Writing out record "<<n<<std::endl;
-      std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-
-      std::stringstream filefn;      filefn << filef << "."<< n;
-      ScidacWriter _ScidacWriter(FGrid->IsBoss());
-      _ScidacWriter.open(filefn.str());
-      _ScidacWriter.writeScidacFieldRecord(src[n],record);
-      _ScidacWriter.close();
-    }
-      
-    FGrid->Barrier();
-
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-    std::cout << GridLogMessage << " Reading back in the single process view "<<std::endl;
-    std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-      
-    for(int n=0;n<nrhs;n++){
-      if ( n==me ) { 
-	std::stringstream filefn;	filefn << filef << "."<< n;
-	ScidacReader  _ScidacReader;
-	_ScidacReader.open(filefn.str());
-	_ScidacReader.readScidacFieldRecord(s_src,record);
-	_ScidacReader.close();
-      }
-    }
     FGrid->Barrier();
   }
 

tests/solver/Test_dwf_mrhs_cg_mpi.cc

@@ -38,7 +38,7 @@ int main (int argc, char ** argv)
   typedef typename DomainWallFermionR::ComplexField ComplexField; 
   typename DomainWallFermionR::ImplParams params; 
 
-  const int Ls=4;
+  const int Ls=8;
 
   Grid_init(&argc,&argv);
 
@@ -69,6 +69,8 @@ int main (int argc, char ** argv)
     }
   }
 
+ 
+  double stp = 1.e-5;
   int nrhs = 1;
   int me;
   for(int i=0;i<mpi_layout.size();i++) nrhs *= (mpi_layout[i]/mpi_split[i]);
@@ -90,7 +92,7 @@ int main (int argc, char ** argv)
   ///////////////////////////////////////////////
   std::vector<int> seeds({1,2,3,4});
 
-  std::vector<FermionField>    src(nrhs,FGrid);
+  std::vector<FermionField> src(nrhs,FGrid);
   std::vector<FermionField> src_chk(nrhs,FGrid);
   std::vector<FermionField> result(nrhs,FGrid);
   FermionField tmp(FGrid);
@@ -123,25 +125,34 @@ int main (int argc, char ** argv)
   for(int s=0;s<nrhs;s++) {
     random(pRNG5,src[s]);
     tmp = 10.0*s;
-    src[s] = (src[s] * 0.1) + tmp;
+//    src[s] = (src[s] * 0.1) + tmp;
     std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
   }
 #endif
   std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
 
   LatticeGaugeField Umu(UGrid); 
-  if(1) { 
+  FieldMetaData header;
+    std::string file("./lat.in");
+    SU3::ColdConfiguration(Umu);
+    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
+  if(0) { 
+    NerscIO::readConfiguration(Umu,header,file);
+    std::cout << GridLogMessage << " "<<file<<" successfully read" <<std::endl;
+  } else {
     GridParallelRNG pRNG(UGrid );  
     std::cout << GridLogMessage << "Intialising 4D RNG "<<std::endl;
     pRNG.SeedFixedIntegers(seeds);
     std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
     SU3::HotConfiguration(pRNG,Umu);
     std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
-    //    std::cout << " Site zero "<< Umu._odata[0]   <<std::endl;
-  } else { 
-    SU3::ColdConfiguration(Umu);
-    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
-  }
+    std::cout << " Site zero "<< Umu._odata[0]   <<std::endl;
+  } 
+   int precision32 = 0;
+   int tworow      = 0;
+   std::string file2("./lat.out");
+   NerscIO::writeConfiguration(Umu,file2,tworow,precision32);
+   std::cout << GridLogMessage << " Successfully saved to " <<file2 <<std::endl;
   /////////////////
   // MPI only sends
   /////////////////
@@ -197,9 +208,9 @@ int main (int argc, char ** argv)
 
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOpCk(Dchk);
-  ConjugateGradient<FermionField> CG((1.0e-2),10000);
+  ConjugateGradient<FermionField> CG((stp),10000);
   s_res = zero;
-  CG(HermOp,s_src,s_res);
+//  CG(HermOp,s_src,s_res);
 
   std::cout << GridLogMessage << " split residual norm "<<norm2(s_res)<<std::endl;
   /////////////////////////////////////////////////////////////
@@ -227,5 +238,15 @@ int main (int argc, char ** argv)
     std::cout << GridLogMessage<<" resid["<<n<<"]  "<< norm2(tmp)/norm2(src[n])<<std::endl;
   }
 
+  for(int s=0;s<nrhs;s++) result[s]=zero;
+
+  int blockDim = 0;//not used for BlockCGVec
+  BlockConjugateGradient<FermionField>    BCGV  (BlockCGVec,blockDim,stp,10000);
+  BCGV.PrintInterval=10;
+{
+  BCGV(HermOpCk,src,result);
+}
+
+
   Grid_finalize();
 }

tests/solver/Test_mobius_bcg.cc

@@ -0,0 +1,277 @@
+   /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_mrhs_cg.cc
+
+    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 */
+#include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main (int argc, char ** argv)
+{
+  typedef typename MobiusFermionR::FermionField FermionField; 
+  typedef typename MobiusFermionR::ComplexField ComplexField; 
+  typename MobiusFermionR::ImplParams params; 
+
+  const int Ls=12;
+
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  std::vector<int> mpi_split (mpi_layout.size(),1);
+  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());
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * rbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  /////////////////////////////////////////////
+  // Split into 1^4 mpi communicators
+  /////////////////////////////////////////////
+
+  for(int i=0;i<argc;i++){
+    if(std::string(argv[i]) == "--split"){
+      for(int k=0;k<mpi_layout.size();k++){
+	std::stringstream ss; 
+	ss << argv[i+1+k]; 
+	ss >> mpi_split[k];
+      }
+      break;
+    }
+  }
+
+ 
+  double stp = 1.e-5;
+  int nrhs = 1;
+  int me;
+  for(int i=0;i<mpi_layout.size();i++){
+//	split_dim[i] = (mpi_layout[i]/mpi_split[i]);
+	nrhs *= (mpi_layout[i]/mpi_split[i]);
+//	split_coor[i] = FGrid._processor_coor[i]/mpi_split[i];
+  }
+  std::cout << GridLogMessage << "Creating split grids " <<std::endl;
+  GridCartesian         * SGrid = new GridCartesian(GridDefaultLatt(),
+						    GridDefaultSimd(Nd,vComplex::Nsimd()),
+						    mpi_split,
+						    *UGrid,me); 
+  std::cout << GridLogMessage <<"Creating split ferm grids " <<std::endl;
+
+  GridCartesian         * SFGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,SGrid);
+  std::cout << GridLogMessage <<"Creating split rb grids " <<std::endl;
+  GridRedBlackCartesian * SrbGrid  = SpaceTimeGrid::makeFourDimRedBlackGrid(SGrid);
+  std::cout << GridLogMessage <<"Creating split ferm rb grids " <<std::endl;
+  GridRedBlackCartesian * SFrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,SGrid);
+  std::cout << GridLogMessage << "Made the grids"<<std::endl;
+  ///////////////////////////////////////////////
+  // Set up the problem as a 4d spreadout job
+  ///////////////////////////////////////////////
+  std::vector<int> seeds({1,2,3,4});
+
+  std::vector<FermionField> src(nrhs,FGrid);
+  std::vector<FermionField> src_chk(nrhs,FGrid);
+  std::vector<FermionField> result(nrhs,FGrid);
+  FermionField tmp(FGrid);
+  std::cout << GridLogMessage << "Made the Fermion Fields"<<std::endl;
+
+  for(int s=0;s<nrhs;s++) result[s]=zero;
+#undef LEXICO_TEST
+#ifdef LEXICO_TEST
+  {
+    LatticeFermion lex(FGrid);  lex = zero;
+    LatticeFermion ftmp(FGrid);
+    Integer stride =10000;
+    double nrm;
+    LatticeComplex coor(FGrid);
+    for(int d=0;d<5;d++){
+      LatticeCoordinate(coor,d);
+      ftmp = stride;
+      ftmp = ftmp * coor;
+      lex = lex + ftmp;
+      stride=stride/10;
+    }
+    for(int s=0;s<nrhs;s++) {
+      src[s]=lex;
+      ftmp = 1000*1000*s;
+      src[s] = src[s] + ftmp;
+    }    
+  }
+#else
+  GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
+  for(int s=0;s<nrhs;s++) {
+    random(pRNG5,src[s]);
+    tmp = 10.0*s;
+//    src[s] = (src[s] * 0.1) + tmp;
+    std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
+  }
+#endif
+  std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
+
+  LatticeGaugeField Umu(UGrid); 
+  FieldMetaData header;
+    std::string file("./lat.in.32IDfine");
+    SU3::ColdConfiguration(Umu);
+    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
+  if(1) { 
+    NerscIO::readConfiguration(Umu,header,file);
+    std::cout << GridLogMessage << " "<<file<<" successfully read" <<std::endl;
+  } else {
+    GridParallelRNG pRNG(UGrid );  
+    std::cout << GridLogMessage << "Intialising 4D RNG "<<std::endl;
+    pRNG.SeedFixedIntegers(seeds);
+    std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
+    SU3::HotConfiguration(pRNG,Umu);
+    std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
+    std::cout << " Site zero "<< Umu._odata[0]   <<std::endl;
+  } 
+   int precision32 = 0;
+   int tworow      = 0;
+   std::string file2("./lat.out");
+   NerscIO::writeConfiguration(Umu,file2,tworow,precision32);
+   std::cout << GridLogMessage << " Successfully saved to " <<file2 <<std::endl;
+  /////////////////
+  // MPI only sends
+  /////////////////
+  LatticeGaugeField s_Umu(SGrid);
+  FermionField s_src(SFGrid);
+  FermionField s_tmp(SFGrid);
+  FermionField s_res(SFGrid);
+
+  std::cout << GridLogMessage << "Made the split grid fields"<<std::endl;
+  ///////////////////////////////////////////////////////////////
+  // split the source out using MPI instead of I/O
+  ///////////////////////////////////////////////////////////////
+  Grid_split  (Umu,s_Umu);
+  Grid_split  (src,s_src);
+  std::cout << GridLogMessage << " split rank  " <<me << " s_src "<<norm2(s_src)<<std::endl;
+
+#ifdef LEXICO_TEST
+  FermionField s_src_tmp(SFGrid);
+  FermionField s_src_diff(SFGrid);
+  {
+    LatticeFermion lex(SFGrid);  lex = zero;
+    LatticeFermion ftmp(SFGrid);
+    Integer stride =10000;
+    double nrm;
+    LatticeComplex coor(SFGrid);
+    for(int d=0;d<5;d++){
+      LatticeCoordinate(coor,d);
+      ftmp = stride;
+      ftmp = ftmp * coor;
+      lex = lex + ftmp;
+      stride=stride/10;
+    }
+    s_src_tmp=lex;
+    ftmp = 1000*1000*me;
+    s_src_tmp = s_src_tmp + ftmp;
+  }
+  s_src_diff = s_src_tmp - s_src;
+  std::cout << GridLogMessage <<" LEXICO test:  s_src_diff " << norm2(s_src_diff)<<std::endl;
+#endif
+
+  ///////////////////////////////////////////////////////////////
+  // Set up N-solvers as trivially parallel
+  ///////////////////////////////////////////////////////////////
+  std::cout << GridLogMessage << " Building the solvers"<<std::endl;
+  RealD mass=0.00107;
+//  RealD mass=0.01;
+  RealD M5=1.8;
+  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,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;
+  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
+
+  MdagMLinearOperator<MobiusFermionR,FermionField> HermOp(Ddwf);
+  MdagMLinearOperator<MobiusFermionR,FermionField> HermOpCk(Dchk);
+  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;
+  /////////////////////////////////////////////////////////////
+  // Report how long they all took
+  /////////////////////////////////////////////////////////////
+  std::vector<uint32_t> iterations(nrhs,0);
+  iterations[me] = CG.IterationsToComplete;
+
+  for(int n=0;n<nrhs;n++){
+    UGrid->GlobalSum(iterations[n]);
+    std::cout << GridLogMessage<<" Rank "<<n<<" "<< iterations[n]<<" CG iterations"<<std::endl;
+  }
+
+  /////////////////////////////////////////////////////////////
+  // Gather and residual check on the results
+  /////////////////////////////////////////////////////////////
+  std::cout << GridLogMessage<< "Unsplitting the result"<<std::endl;
+  Grid_unsplit(result,s_res);
+
+
+  std::cout << GridLogMessage<< "Checking the residuals"<<std::endl;
+  for(int n=0;n<nrhs;n++){
+    std::cout << GridLogMessage<< " res["<<n<<"] norm "<<norm2(result[n])<<std::endl;
+    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,100000);
+  BCGV.PrintInterval=10;
+{
+  BCGV(HermOpCk,src,result);
+}
+
+
+  Grid_finalize();
+}

tests/solver/Test_pcg.cc

@@ -1,247 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_pcg.cc
-
-    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 */
-#include<bitset>
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-//Preconditioning:   M psi = chi
-//               =  M P^-1 P psi = chi
-//               =  M P^-1 psi' = chi
-
-//Solve for psi' using M P^-1 as operator, then apply  P^-1 psi' = psi 
-
-//Inexact preconditioned CG requires slight modification because we want to avoid computing P^-1 exactly
-
-
-/////////////////////////////////////////////////////////////
-// Base classes for iterative processes based on operators
-// single input vec, single output vec.
-/////////////////////////////////////////////////////////////
-
-template<typename T>
-T parse(const std::string &name, std::istream &in){
-  std::string p;
-  in >> p;
-  assert(p==name);
-  char eq;
-  in >> eq;
-  assert(eq == '=');
-  T out;
-  in >> out;
-  return out;
-}
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  int Ls=12;
-  RealD mass=0.01;
-  RealD outer_tol = 1e-8;
-  RealD inner_tol_full = 1e-5;
-  RealD inner_tol_half = 1e-5;
-  RealD inner_tol_16c = 1e-5;
-  RealD inner_tol_8c = 1e-5;
-
-  RealD relup_delta_full = 0.1;
-  RealD relup_delta_half = 0.1;
-  RealD relup_delta_16c = 0.1;
-  RealD relup_delta_8c = 0.1;
-
-  std::string config_file = "";
-  RealD lo = 1.0;
-  RealD hi = 64.0;
-  int   order=10;
-  for(int i=1;i<argc;i++){
-    if(std::string(argv[i]) == "--params"){
-      std::ifstream f(argv[i+1]);
-      f.exceptions ( std::ifstream::failbit | std::ifstream::badbit );
-      Ls = parse<int>("Ls", f);
-#define PARSEIT(NM) NM = parse<RealD>(#NM, f)
-      PARSEIT(mass);
-      PARSEIT(outer_tol);
-      PARSEIT(inner_tol_full);
-      PARSEIT(inner_tol_half);
-      PARSEIT(inner_tol_16c);
-      PARSEIT(inner_tol_8c);
-      PARSEIT(relup_delta_full);
-      PARSEIT(relup_delta_half);
-      PARSEIT(relup_delta_16c);
-      PARSEIT(relup_delta_8c);
-#undef PARSEIT
-
-      //f >> outer_tol >> inner_tol_full >> inner_tol_half >> inner_tol_16c >> inner_tol_8c;      
-    }else if(std::string(argv[i]) == "--config"){
-      config_file = argv[i+1];
-    }else if(std::string(argv[i]) == "--order"){
-      std::string ss(argv[i+1]);
-      std::stringstream f(ss);
-      f>>order; std::cout << " Order poly set to " <<order<<std::endl;
-    }else if(std::string(argv[i]) == "--lo"){
-      std::string ss(argv[i+1]);
-      std::stringstream f(ss);
-      f>>lo; std::cout << " Lo poly set to " <<lo<<std::endl;
-    }else if(std::string(argv[i]) == "--hi"){
-      std::string ss(argv[i+1]);
-      std::stringstream f(ss);
-      f>>hi; std::cout << " Hi poly set to " <<hi<<std::endl;
-    }
-  }
-  
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  GridCartesian         * UGrid_f   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f);
-  GridCartesian         * FGrid_f   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid_f);
-  GridRedBlackCartesian * FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid_f);
-  
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
-
-  LatticeFermionD    src(FGrid); random(RNG5,src);
-  LatticeFermionD result(FGrid); result=zero;
-  LatticeGaugeFieldD Umu(UGrid);
-  LatticeGaugeFieldF Umu_f(UGrid_f); 
-
-  if(config_file.size() > 0){
-    FieldMetaData header;
-    NerscIO::readConfiguration(Umu,header,config_file);
-  }else{
-    SU3::HotConfiguration(RNG4,Umu);
-  }
-
-  precisionChange(Umu_f,Umu);
-  
-  RealD M5=1.8;
-
-  LatticeFermionD    src_o(FrbGrid);
-  pickCheckerboard(Odd,src_o,src);
-
-  //if(0){ //Test preconditioned CG
-    std::cout << "Test preconditioned CG" << std::endl;
-    LatticeFermionD result_o(FrbGrid);
-    LatticeFermionD result_o_2(FrbGrid);
-    result_o.checkerboard = Odd;
-    result_o = zero;
-    result_o_2.checkerboard = Odd;
-    result_o_2 = zero;
-
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-    //DoNothingLinearOperator<LatticeFermionD> Prec;
-    //FixedIterConjugateGradientPreconditioner<LatticeFermionD> Prec(HermOpEO, 20);
-    //    SloppyConjugateGradientPreconditioner<LatticeFermionD> Prec(HermOpEO, 1e-2, 1000);
-    PolynomialPreconditioner<LatticeFermionD> Prec(HermOpEO,lo,hi,order) ;
-
-    std::cout << "Preconditioned CG" << std::endl;
-    InexactPreconditionedConjugateGradient<LatticeFermionD> pCG(Prec,1.0e-8,10000);
-    pCG(HermOpEO,src_o,result_o);
-
-    std::cout << "Starting regular CG" << std::endl;
-    ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-    CG(HermOpEO,src_o,result_o_2);
-
-    LatticeFermionD diff_o(FrbGrid);
-    RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
-
-    std::cout << "pCG HermOp applications " << " Lo " << lo << " Hi " << hi << " Order " << order << " " <<  pCG.IterationsToComplete << "(outer) + " << Prec.InnerIterations << "(inner) = " << pCG.IterationsToComplete + Prec.InnerIterations << std::endl;
-    std::cout << "CG HermOp applications " << CG.IterationsToComplete << std::endl;
-    std::cout << "Diff between results: " << diff << std::endl;
-  //}
-
-  if(0){ //Test compressor
-    LatticeFermionD result_o(FrbGrid);
-    LatticeFermionD result_o_2(FrbGrid);
-    result_o.checkerboard = Odd;
-    result_o = zero;
-    result_o_2.checkerboard = Odd;
-    result_o_2 = zero;
-
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-
-    DomainWallFermionDF DdwfC(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionDF,LatticeFermionD> HermOpEOC(DdwfC);
-
-    std::cout << "Starting regular CG with compressed operator" << std::endl;
-    Integer iter1;
-    {
-      ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-      CG.ErrorOnNoConverge = false;
-      CG(HermOpEOC,src_o,result_o);
-      iter1 = CG.IterationsToComplete;
-    }
-    Integer iter2;
-    {
-      std::cout << "Starting regular CG" << std::endl;
-      ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-      CG(HermOpEO,src_o,result_o_2);
-      iter2 = CG.IterationsToComplete;
-    }
-
-    LatticeFermionD diff_o(FrbGrid);
-    RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
-
-    std::cout << "CG HermOp CC applications " << iter1 << std::endl;
-    std::cout << "CG HermOp applications " << iter2 << std::endl;
-    std::cout << "Diff between results: " << diff << std::endl;
-  }
-  
-  if(1){ //Compare mixed prec restarted single/single internal with same but with single/compressed
-    LatticeFermionD result_o_full(FrbGrid);
-    LatticeFermionD result_o_half(FrbGrid);
-    LatticeFermionD result_o_16(FrbGrid);
-    LatticeFermionD result_o_8(FrbGrid);
-    result_o_full.checkerboard = Odd;
-    result_o_full = zero;
-    result_o_16 = result_o_8 = result_o_half = result_o_full;
-
-    //Std
-    DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-
-    DomainWallFermionF Ddwf_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionF,LatticeFermionF> HermOpEO_f(Ddwf_f);
-
-    //1/2 prec
-    DomainWallFermionFH Ddwfhalf_f(Umu_f,*FGrid_f,*FrbGrid_f,*UGrid_f,*UrbGrid_f,mass,M5);
-    SchurDiagMooeeOperator<DomainWallFermionFH,LatticeFermionF> HermOpEOhalf_f(Ddwfhalf_f);
-
-  }
-
-
-  Grid_finalize();
-}