diff --git a/Grid/algorithms/blas/BatchedBlas.h b/Grid/algorithms/blas/BatchedBlas.h
index 3f4856f0..e376bb18 100644
--- a/Grid/algorithms/blas/BatchedBlas.h
+++ b/Grid/algorithms/blas/BatchedBlas.h
@@ -270,7 +270,6 @@ public:
     assert(err==CUBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_SYCL
-    //    std::cerr << " Calling SYCL batched ZGEMM "<<std::endl;
       int64_t m64=m;
       int64_t n64=n;
       int64_t k64=k;
@@ -506,7 +505,7 @@ public:
 						  (ComplexF *) &beta_p[0],
 						  (ComplexF **)&Cmn[0], (const int64_t *)&ldc64,
 						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
-      synchronise();
+    synchronise();
 #endif
 #if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
     // Need a default/reference implementation; use Eigen
diff --git a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h
index 8168aa9d..66d2812d 100644
--- a/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h
+++ b/Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h
@@ -279,11 +279,11 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
       diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
       _sort.push(eval2,Nm);
-      //      Glog << "#Ritz value before shift: "<< std::endl;
+      Glog << "#Ritz value before shift: "<< std::endl;
       for(int i=0; i<Nm; ++i){
-	//        std::cout.precision(13);
-	//        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-	//        std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+	std::cout.precision(13);
+	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
       }
       
       //----------------------------------------------------------------------
@@ -297,7 +297,8 @@ public:
         
         unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
 
-        for(int ip=Nk; ip<Nm; ++ip){ 
+        for(int ip=Nk; ip<Nm; ++ip){
+	  Glog << " ip "<<ip<<" / "<<Nm<<std::endl;
           shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
         }
         
@@ -325,7 +326,7 @@ public:
         Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
         diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
         _sort.push(eval2,Nk);
-	//        Glog << "#Ritz value after shift: "<< std::endl;
+	Glog << "#Ritz value after shift: "<< std::endl;
         for(int i=0; i<Nk; ++i){
 	  //          std::cout.precision(13);
 	  //          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
@@ -467,10 +468,10 @@ public:
   
     // set initial vector
     for (int i=0; i<Nu; ++i) {
-      //      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
       evec[i] = src[i];
       orthogonalize(evec[i],evec,i);
-      //      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
+      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
     }
 //    exit(-43);
     
@@ -506,11 +507,11 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nr,Nr);
       diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid);
       _sort.push(eval2,Nr);
-      //      Glog << "#Ritz value: "<< std::endl;
+      Glog << "#Ritz value: "<< std::endl;
       for(int i=0; i<Nr; ++i){
-	//        std::cout.precision(13);
-	//        std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-	//        std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+	std::cout.precision(13);
+	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
       }
       
       // Convergence test
@@ -570,6 +571,7 @@ public:
       Glog << fname + " NOT converged ; Summary :\n";
     } else {
       Glog << fname + " CONVERGED ; Summary :\n";
+      Nstop = Nconv_guess; // Just take them all
       // Sort convered eigenpairs.
       std::vector<Field>  Btmp(Nstop,grid); // waste of space replicating
 
@@ -642,7 +644,7 @@ private:
       //      for (int u=0; u<mrhs; ++u) Glog << " out["<<u<<"] = "<<norm2(out[u])<<std::endl;
       k_start +=mrhs;
     }
-    //    Glog << "LinAlg "<< std::endl;
+    Glog << "LinAlg "<< std::endl;
     
     if (b>0) {
       for (int u=0; u<Nu; ++u) {
@@ -676,7 +678,7 @@ private:
       }
       w_copy[u] = w[u];
     }
-    //    Glog << "LinAlg done"<< std::endl;
+    Glog << "LinAlg done"<< std::endl;
     
     // In block version, the steps 6 and 7 in Lanczos construction is
     // replaced by the QR decomposition of new basis block.
@@ -689,15 +691,15 @@ private:
     }
 
     // re-orthogonalization for numerical stability
-    //    Glog << "Gram Schmidt"<< std::endl;
+    Glog << "Gram Schmidt"<< std::endl;
     orthogonalize(w,Nu,evec,R);
     // QR part
     for (int u=1; u<Nu; ++u) {
       orthogonalize(w[u],w,u);
     }
-    //    Glog << "Gram Schmidt done "<< std::endl;
+    Glog << "Gram Schmidt done "<< std::endl;
     
-    //    Glog << "LinAlg "<< std::endl;
+    Glog << "LinAlg "<< std::endl;
     for (int u=0; u<Nu; ++u) {
       //for (int v=0; v<Nu; ++v) {
       for (int v=u; v<Nu; ++v) {
@@ -714,7 +716,7 @@ private:
 	//        Glog <<" In block "<< b << "," <<" beta[" << u << "," << k-L << "] = " << lme[u][k] << std::endl;
       }
     }
-    //    Glog << "LinAlg done "<< std::endl;
+    Glog << "LinAlg done "<< std::endl;
 
     if (b < Nm/Nu-1) {
       for (int u=0; u<Nu; ++u) {
@@ -779,7 +781,7 @@ private:
     
     for ( int u=0; u<Nu; ++u ) {
       for (int k=0; k<Nk; ++k ) {
-//        Glog << "lmd "<<u<<" "<<k<<" "<<lmd[u][k] -conjugate(lmd[u][k])<<std::endl;
+	//	Glog << "lmd "<<u<<" "<<k<<" "<<lmd[u][k] -conjugate(lmd[u][k])<<std::endl;
         BlockTriDiag(k,u+(k/Nu)*Nu) = lmd[u][k];
       }
     }
@@ -933,7 +935,7 @@ if (1){
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    //Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
+    Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     M = Eigen::MatrixXcd::Zero(Nk,Nk);
@@ -951,7 +953,7 @@ if (1){
         M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
       }
     }
-    //Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
+    Glog << "unpackHermitBlockTriDiagMatToEigen() end" << std::endl; 
   }
  
 
@@ -961,7 +963,7 @@ if (1){
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    //Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
+    Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     
@@ -977,7 +979,7 @@ if (1){
         lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
       }
     }
-    //Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
+    Glog << "packHermitBlockTriDiagMatfromEigen() end" <<std::endl; 
   }
 
 
@@ -986,7 +988,7 @@ if (1){
 		            RealD Dsh,
 		            Eigen::MatrixXcd& Qprod)
   {
-    //Glog << "shiftedQRDecompEigen() begin" << '\n'; 
+    Glog << "shiftedQRDecompEigen() begin" << '\n'; 
     Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
     Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm);
     Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
@@ -1002,6 +1004,7 @@ if (1){
                         // lower triangular part used to represent series
                         // of Q sequence.
 
+    Glog << "shiftedQRDecompEigen() Housholder & QR" << '\n'; 
     // equivalent operation of Qprod *= Q
     //M = Eigen::MatrixXcd::Zero(Nm,Nm);
     
@@ -1022,6 +1025,7 @@ if (1){
     
     Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
 
+    Glog << "shiftedQRDecompEigen() Mtmp create" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=0; j<Nm-(Nu+1); ++j) {
         for (int k=0; k<Nu+1+j; ++k) {
@@ -1029,6 +1033,7 @@ if (1){
         }
       }
     }
+    Glog << "shiftedQRDecompEigen() Mtmp loop1" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=Nm-(Nu+1); j<Nm; ++j) {
         for (int k=0; k<Nm; ++k) {
@@ -1036,6 +1041,7 @@ if (1){
         }
       }
     }
+    Glog << "shiftedQRDecompEigen() Mtmp loop2" << '\n'; 
     
     //static int ntimes = 2;
     //for (int j=0; j<Nm-(ntimes*Nu); ++j) {
@@ -1061,11 +1067,13 @@ if (1){
         Mtmp(j,i) = conj(Mtmp(i,j));
       }
     }
+    Glog << "shiftedQRDecompEigen() Mtmp loop3" << '\n'; 
 
     for (int i=0; i<Nm; ++i) {
       Mtmp(i,i) = real(Mtmp(i,i)) + Dsh;
     }
     
+    Glog << "shiftedQRDecompEigen() Mtmp loop4" << '\n'; 
     M = Mtmp;
 
     //M = Q.adjoint()*(M*Q);
@@ -1077,7 +1085,7 @@ if (1){
     //  }
     //}
     
-    //Glog << "shiftedQRDecompEigen() end" << endl; 
+    Glog << "shiftedQRDecompEigen() end" <<std::endl; 
   }
 
   void exampleQRDecompEigen(void)
diff --git a/Grid/algorithms/iterative/SchurRedBlack.h b/Grid/algorithms/iterative/SchurRedBlack.h
index d97e4993..494aa77b 100644
--- a/Grid/algorithms/iterative/SchurRedBlack.h
+++ b/Grid/algorithms/iterative/SchurRedBlack.h
@@ -499,6 +499,87 @@ namespace Grid {
       }
   };
 
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Site diagonal is identity, left preconditioned by Mee^inv
+  // ( 1 - Mee^inv Meo Moo^inv Moe ) phi = Mee_inv ( Mee - Meo Moo^inv Moe Mee^inv  ) phi =  Mee_inv eta
+  //
+  // Solve:
+  // ( 1 - Mee^inv Meo Moo^inv Moe )^dag ( 1 - Mee^inv Meo Moo^inv Moe ) phi = ( 1 - Mee^inv Meo Moo^inv Moe )^dag  Mee_inv eta
+  //
+  // Old notation e<->o
+  //
+  // Left precon by Moo^-1
+  //  b) (Doo^{dag} M_oo^-dag) (Moo^-1 Doo) psi_o =  [ (D_oo)^dag M_oo^-dag ] Moo^-1 L^{-1}  eta_o
+  //                                   eta_o'     = (D_oo)^dag  M_oo^-dag Moo^-1 (eta_o - Moe Mee^{-1} eta_e)
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  template<class Field> class SchurRedBlackDiagOneSolve : public SchurRedBlackBase<Field> {
+  public:
+    typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
+
+    /////////////////////////////////////////////////////
+    // Wrap the usual normal equations Schur trick
+    /////////////////////////////////////////////////////
+  SchurRedBlackDiagOneSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false,
+      const bool _solnAsInitGuess = false)  
+    : SchurRedBlackBase<Field>(HermitianRBSolver,initSubGuess,_solnAsInitGuess) {};
+
+    virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)
+    {
+      GridBase *grid = _Matrix.RedBlackGrid();
+      GridBase *fgrid= _Matrix.Grid();
+
+      SchurDiagOneOperator<Matrix,Field> _HermOpEO(_Matrix);
+      
+      Field   tmp(grid);
+      Field  Mtmp(grid);
+
+      pickCheckerboard(Even,src_e,src);
+      pickCheckerboard(Odd ,src_o,src);
+    
+      /////////////////////////////////////////////////////
+      // src_o = Mpcdag *MooeeInv * (source_o - Moe MeeInv source_e)
+      /////////////////////////////////////////////////////
+      _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.Checkerboard() ==Even);
+      _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.Checkerboard() ==Odd);     
+      Mtmp=src_o-Mtmp;                 
+      _Matrix.MooeeInv(Mtmp,tmp);      assert( tmp.Checkerboard() ==Odd);     
+      
+      // get the right MpcDag
+      _HermOpEO.MpcDag(tmp,src_o);     assert(src_o.Checkerboard() ==Odd);       
+    }
+
+    virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol)
+    {
+      GridBase *grid = _Matrix.RedBlackGrid();
+      GridBase *fgrid= _Matrix.Grid();
+
+      Field   tmp(grid);
+      Field   sol_e(grid);
+
+
+      ///////////////////////////////////////////////////
+      // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+      ///////////////////////////////////////////////////
+      _Matrix.Meooe(sol_o,tmp);    assert(  tmp.Checkerboard()   ==Even);
+      tmp = src_e-tmp;             assert(  src_e.Checkerboard() ==Even);
+      _Matrix.MooeeInv(tmp,sol_e); assert(  sol_e.Checkerboard() ==Even);
+     
+      setCheckerboard(sol,sol_e);  assert(  sol_e.Checkerboard() ==Even);
+      setCheckerboard(sol,sol_o);  assert(  sol_o.Checkerboard() ==Odd );
+    };
+
+    virtual void RedBlackSolve   (Matrix & _Matrix,const Field &src_o, Field &sol_o)
+    {
+      SchurDiagOneOperator<Matrix,Field> _HermOpEO(_Matrix);
+      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);
+    };
+    virtual void RedBlackSolve   (Matrix & _Matrix,const std::vector<Field> &src_o,  std::vector<Field> &sol_o)
+    {
+      SchurDiagOneOperator<Matrix,Field> _HermOpEO(_Matrix);
+      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); 
+    }
+  };
+
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
   // Site diagonal is identity, right preconditioned by Mee^inv
   // ( 1 - Meo Moo^inv Moe Mee^inv  ) phi =( 1 - Meo Moo^inv Moe Mee^inv  ) Mee psi =  = eta  = eta
diff --git a/Grid/allocator/AlignedAllocator.h b/Grid/allocator/AlignedAllocator.h
index 8a27f527..293ce2fb 100644
--- a/Grid/allocator/AlignedAllocator.h
+++ b/Grid/allocator/AlignedAllocator.h
@@ -54,6 +54,9 @@ public:
     size_type bytes = __n*sizeof(_Tp);
     profilerAllocate(bytes);
     _Tp *ptr = (_Tp*) MemoryManager::CpuAllocate(bytes);
+    if ( (_Tp*)ptr == (_Tp *) NULL ) {
+      printf("Grid CPU Allocator got NULL for %lu bytes\n",(unsigned long) bytes );
+    }
     assert( ( (_Tp*)ptr != (_Tp *)NULL ) );
     return ptr;
   }
@@ -100,6 +103,9 @@ public:
     size_type bytes = __n*sizeof(_Tp);
     profilerAllocate(bytes);
     _Tp *ptr = (_Tp*) MemoryManager::SharedAllocate(bytes);
+    if ( (_Tp*)ptr == (_Tp *) NULL ) {
+      printf("Grid Shared Allocator got NULL for %lu bytes\n",(unsigned long) bytes );
+    }
     assert( ( (_Tp*)ptr != (_Tp *)NULL ) );
     return ptr;
   }
@@ -145,6 +151,9 @@ public:
     size_type bytes = __n*sizeof(_Tp);
     profilerAllocate(bytes);
     _Tp *ptr = (_Tp*) MemoryManager::AcceleratorAllocate(bytes);
+    if ( (_Tp*)ptr == (_Tp *) NULL ) {
+      printf("Grid Device Allocator got NULL for %lu bytes\n",(unsigned long) bytes );
+    }
     assert( ( (_Tp*)ptr != (_Tp *)NULL ) );
     return ptr;
   }
diff --git a/Grid/allocator/MemoryManager.h b/Grid/allocator/MemoryManager.h
index af3b2c8b..e2d7b49e 100644
--- a/Grid/allocator/MemoryManager.h
+++ b/Grid/allocator/MemoryManager.h
@@ -212,6 +212,7 @@ private:
 #endif
 
  public:
+  static void DisplayMallinfo(void);
   static void NotifyDeletion(void * CpuPtr);
   static void Print(void);
   static void PrintAll(void);
diff --git a/Grid/lattice/Lattice_reduction.h b/Grid/lattice/Lattice_reduction.h
index a0abddd4..af62cefe 100644
--- a/Grid/lattice/Lattice_reduction.h
+++ b/Grid/lattice/Lattice_reduction.h
@@ -356,7 +356,8 @@ axpby_norm_fast(Lattice<vobj> &z,sobj a,sobj b,const Lattice<vobj> &x,const Latt
   nrm = real(TensorRemove(sum(inner_tmp_v,sites)));
 #else
   typedef decltype(innerProduct(x_v[0],y_v[0])) inner_t;
-  Vector<inner_t> inner_tmp(sites);
+  deviceVector<inner_t> inner_tmp;
+  inner_tmp.resize(sites);
   auto inner_tmp_v = &inner_tmp[0];
 
   accelerator_for( ss, sites, nsimd,{
diff --git a/Grid/qcd/hmc/GenericHMCrunner.h b/Grid/qcd/hmc/GenericHMCrunner.h
index 1429d848..b53755aa 100644
--- a/Grid/qcd/hmc/GenericHMCrunner.h
+++ b/Grid/qcd/hmc/GenericHMCrunner.h
@@ -90,6 +90,7 @@ public:
         exit(1);
       }
       Parameters.StartingType = arg;
+      std::cout <<GridLogMessage << " GenericHMCrunner --StartingType "<<arg<<std::endl;
     }
 
     if (GridCmdOptionExists(argv, argv + argc, "--StartingTrajectory")) {
@@ -97,6 +98,7 @@ public:
       std::vector<int> ivec(0);
       GridCmdOptionIntVector(arg, ivec);
       Parameters.StartTrajectory = ivec[0];
+      std::cout <<GridLogMessage << " GenericHMCrunner --StartingTrajectory "<<ivec[0]<<std::endl;
     }
 
     if (GridCmdOptionExists(argv, argv + argc, "--Trajectories")) {
@@ -104,6 +106,7 @@ public:
       std::vector<int> ivec(0);
       GridCmdOptionIntVector(arg, ivec);
       Parameters.Trajectories = ivec[0];
+      std::cout << GridLogMessage<<" GenericHMCrunner Command Line --Trajectories "<<ivec[0]<<std::endl;
     }
 
     if (GridCmdOptionExists(argv, argv + argc, "--Thermalizations")) {
@@ -111,6 +114,7 @@ public:
       std::vector<int> ivec(0);
       GridCmdOptionIntVector(arg, ivec);
       Parameters.NoMetropolisUntil = ivec[0];
+      std::cout << GridLogMessage<<" GenericHMCrunner --Thermalizations "<<ivec[0]<<std::endl;
     }
     if (GridCmdOptionExists(argv, argv + argc, "--ParameterFile")) {
       arg = GridCmdOptionPayload(argv, argv + argc, "--ParameterFile");
diff --git a/Grid/qcd/hmc/integrators/Integrator.h b/Grid/qcd/hmc/integrators/Integrator.h
index 0276b6fd..549920a0 100644
--- a/Grid/qcd/hmc/integrators/Integrator.h
+++ b/Grid/qcd/hmc/integrators/Integrator.h
@@ -137,9 +137,11 @@ public:
 
       double start_force = usecond();
 
+      MemoryManager::Print();
       as[level].actions.at(a)->deriv_timer_start();
       as[level].actions.at(a)->deriv(Smearer, force);  // deriv should NOT include Ta
       as[level].actions.at(a)->deriv_timer_stop();
+      MemoryManager::Print();
 
       auto name = as[level].actions.at(a)->action_name();
 
@@ -246,7 +248,11 @@ public:
     }
   };
 
-  virtual ~Integrator() {}
+  virtual ~Integrator()
+  {
+    // Pain in the ass to clean up the Level pointers
+    // Guido's design is at fault as per comment above in constructor
+  }
 
   virtual std::string integrator_name() = 0;
   
@@ -460,6 +466,7 @@ public:
     for (int level = 0; level < as.size(); ++level) {
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
 
+	MemoryManager::Print();
         // get gauge field from the SmearingPolicy and
         // based on the boolean is_smeared in actionID
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] action eval " << std::endl;
@@ -468,6 +475,7 @@ public:
    	        as[level].actions.at(actionID)->S_timer_stop();
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] H = " << Hterm << std::endl;
         H += Hterm;
+	MemoryManager::Print();
 
       }
       as[level].apply(S_hireps, Representations, level, H);
diff --git a/HMC/Mobius2p1f.cc b/HMC/Mobius2p1f.cc
index 4ab1f20f..8042d6e6 100644
--- a/HMC/Mobius2p1f.cc
+++ b/HMC/Mobius2p1f.cc
@@ -58,7 +58,7 @@ int main(int argc, char **argv) {
   HMCparameters HMCparams;
   HMCparams.StartTrajectory  = 0;
   HMCparams.Trajectories     = 200;
-  HMCparams.NoMetropolisUntil=  20;
+  HMCparams.NoMetropolisUntil=  0;
   // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
   HMCparams.StartingType     =std::string("ColdStart");
   HMCparams.MD = MD;
@@ -70,7 +70,7 @@ int main(int argc, char **argv) {
   CheckpointerParameters CPparams;
   CPparams.config_prefix = "ckpoint_EODWF_lat";
   CPparams.rng_prefix    = "ckpoint_EODWF_rng";
-  CPparams.saveInterval  = 10;
+  CPparams.saveInterval  = 1;
   CPparams.format        = "IEEE64BIG";
   TheHMC.Resources.LoadNerscCheckpointer(CPparams);
 
@@ -186,6 +186,8 @@ int main(int argc, char **argv) {
 
   /////////////////////////////////////////////////////////////
   // HMC parameters are serialisable
+  TheHMC.ReadCommandLine(argc,argv);  // params on CML or from param file                                                                     
+  TheHMC.initializeGaugeFieldAndRNGs(U);
 
   std::cout << GridLogMessage << " Running the HMC "<< std::endl;
   TheHMC.Run();  // no smearing
diff --git a/systems/Aurora/sourceme.sh b/systems/Aurora/sourceme.sh
index 98f618ab..7abe667f 100644
--- a/systems/Aurora/sourceme.sh
+++ b/systems/Aurora/sourceme.sh
@@ -24,13 +24,12 @@ export INTELGT_AUTO_ATTACH_DISABLE=1
 # -ftarget-register-alloc-mode=pvc:small
 # -ftarget-register-alloc-mode=pvc:large
 # -ftarget-register-alloc-mode=pvc:auto
-#
+#export MPIR_CVAR_CH4_OFI_ENABLE_HMEM=1
 
 export HTTP_PROXY=http://proxy.alcf.anl.gov:3128
 export HTTPS_PROXY=http://proxy.alcf.anl.gov:3128
 export http_proxy=http://proxy.alcf.anl.gov:3128
 export https_proxy=http://proxy.alcf.anl.gov:3128
-#export MPIR_CVAR_CH4_OFI_ENABLE_HMEM=1
 git config --global http.proxy http://proxy.alcf.anl.gov:3128
 
 #source ~/spack/share/spack/setup-env.sh
diff --git a/systems/Frontier/sourceme.sh b/systems/Frontier/sourceme.sh
index 652e5b45..a6f49d8c 100644
--- a/systems/Frontier/sourceme.sh
+++ b/systems/Frontier/sourceme.sh
@@ -3,7 +3,7 @@ spack load c-lime
 module load emacs 
 module load PrgEnv-gnu
 module load rocm
-module load cray-mpich/8.1.23
+module load cray-mpich
 module load gmp
 module load cray-fftw
 module load craype-accel-amd-gfx90a
diff --git a/tests/debug/Test_8888.cc b/tests/debug/Test_8888.cc
new file mode 100644
index 00000000..279bc8bd
--- /dev/null
+++ b/tests/debug/Test_8888.cc
@@ -0,0 +1,118 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_general_coarse_hdcg.cc
+
+    Copyright (C) 2023
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    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/ImplicitlyRestartedBlockLanczos.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
+#include <Grid/algorithms/iterative/AdefMrhs.h>
+
+using namespace std;
+using namespace Grid;
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls=8;
+  const int nbasis = 40;
+  const int cb = 0 ;
+  RealD mass=0.01;
+  RealD M5=1.8;
+  RealD b=1.0;
+  RealD c=0.0;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  ///////////////////////// RNGs /////////////////////////////////
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  std::vector<int> cseeds({5,6,7,8});
+
+  GridParallelRNG          RNG5(FGrid);   RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
+
+  ///////////////////////// Configuration /////////////////////////////////
+  LatticeGaugeField Umu(UGrid);
+
+  FieldMetaData header;
+  std::string file("ckpoint_EODWF_lat.125");
+  NerscIO::readConfiguration(Umu,header,file);
+
+  //////////////////////// Fermion action //////////////////////////////////
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+
+  MdagMLinearOperator<MobiusFermionD, LatticeFermion> HermOp(Ddwf);
+
+  
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "         Fine Power method            "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+
+  LatticeFermionD pm_src(FGrid);
+  pm_src = ComplexD(1.0);
+  PowerMethod<LatticeFermionD>       fPM;
+  fPM(HermOp,pm_src);
+
+  
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "         Fine Lanczos  (poly, low)    "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  
+  int Nk=80;
+  int Nm=Nk*3;
+  int Nstop=8;
+  int Nconv_test_interval=1;
+  
+  //  Chebyshev<LatticeFermionD>      IRLChebyLo(0.2,64.0,201);  // 1 iter
+  Chebyshev<LatticeFermionD>      IRLChebyLo(0.0,55.0,101);  // 1 iter
+  FunctionHermOp<LatticeFermionD>    PolyOp(IRLChebyLo,HermOp);
+  PlainHermOp<LatticeFermionD>          Op(HermOp);
+
+  ImplicitlyRestartedLanczos IRL(PolyOp,
+				 Op,
+				 Nk, // sought vecs
+				 Nk, // sought vecs
+				 Nm, // spare vecs
+				 1.0e-8,
+				 10 // Max iterations
+				 );
+
+  int Nconv;
+  std::vector<RealD>            eval(Nm);
+  std::vector<LatticeFermionD>     evec(Nm,FGrid);
+  LatticeFermionD     irl_src(FGrid);
+
+  IRL.calc(eval,evec,irl_src,Nconv);
+
+  Grid_finalize();
+  return 0;
+}
diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc
index 74492fd9..068c260f 100644
--- a/tests/debug/Test_cayley_cg.cc
+++ b/tests/debug/Test_cayley_cg.cc
@@ -392,9 +392,27 @@ void  TestCGschur(What & Ddwf,
 		   GridParallelRNG *RNG5)
 {
   LatticeFermion src   (FGrid); random(*RNG5,src);
-  LatticeFermion result(FGrid); result=Zero();
+  LatticeFermion result1(FGrid); result1=Zero();
+  LatticeFermion result2(FGrid); result2=Zero();
+  LatticeFermion result3(FGrid); result3=Zero();
 
   ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
   SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
-  SchurSolver(Ddwf,src,result);
+  SchurSolver(Ddwf,src,result1);
+
+  SchurRedBlackDiagOneSolve<LatticeFermion> SchurSolverSymm1(CG);
+  SchurSolverSymm1(Ddwf,src,result2);
+
+  SchurRedBlackDiagTwoSolve<LatticeFermion> SchurSolverSymm2(CG);
+  SchurSolverSymm2(Ddwf,src,result3);
+
+  std::cout << GridLogMessage << " Standard " <<norm2(result1)<<std::endl;
+
+  std::cout << GridLogMessage << " Symm1    " <<norm2(result2)<<std::endl; 
+  result2=result2-result1;
+  std::cout << GridLogMessage << " diff " <<norm2(result2) <<std::endl; 
+
+  std::cout << GridLogMessage << " Symm2    " <<norm2(result3)<<std::endl; 
+  result3=result3-result1;
+  std::cout << GridLogMessage << " diff " <<norm2(result3) <<std::endl; 
 }
diff --git a/tests/debug/Test_general_coarse.cc b/tests/debug/Test_general_coarse.cc
index 28130a3c..4351a901 100644
--- a/tests/debug/Test_general_coarse.cc
+++ b/tests/debug/Test_general_coarse.cc
@@ -244,7 +244,7 @@ int main (int argc, char ** argv)
 
   GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); 
 
-  
+#if 0  
   MultiGeneralCoarsenedMatrix mrhs(LittleDiracOp,CoarseMrhs);
   typedef decltype(mrhs) MultiGeneralCoarsenedMatrix_t;
   
@@ -307,7 +307,8 @@ int main (int argc, char ** argv)
     rh_res= Zero();
     mrhsCG(MrhsCoarseOp,rh_src,rh_res);
   }
-  
+
+#endif
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
diff --git a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc
index c464b25b..054b4c66 100644
--- a/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc
+++ b/tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc
@@ -145,7 +145,7 @@ int main (int argc, char ** argv)
   Grid_init(&argc,&argv);
 
   const int Ls=24;
-  const int nbasis = 60;
+  const int nbasis = 62;
   const int cb = 0 ;
   RealD mass=0.00078;
   RealD M5=1.8;
@@ -160,7 +160,7 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
   // Construct a coarsened grid with 4^4 cell
-  Coordinate Block({4,4,4,4});
+  Coordinate Block({4,4,6,4});
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/Block[d];
diff --git a/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc b/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc
index 71ad42ba..c45b2cb1 100644
--- a/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc
+++ b/tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc
@@ -160,7 +160,8 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
   // Construct a coarsened grid with 4^4 cell
-  Coordinate Block({4,4,6,6});
+  //  Coordinate Block({4,4,6,4});
+  Coordinate Block({4,4,4,4});
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/Block[d];
@@ -217,7 +218,7 @@ int main (int argc, char ** argv)
   std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/evecs.scidac");
   std::string eval_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/eval.xml");
   bool load_agg=true;
-  bool load_refine=false;
+  bool load_refine=true;
   bool load_mat=false;
   bool load_evec=false;
 
@@ -276,17 +277,25 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
 
   typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
-  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,301);  // 1 iter
+  //  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,301);  // 4.4.6.4
+  //  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,501);  // for 4.4.4.4 blocking 350 evs
+  //  Chebyshev<CoarseVector>      IRLCheby(0.0014,42.0,501);  // for 4.4.4.4 blocking 700 evs
+  //  Chebyshev<CoarseVector>      IRLCheby(0.002,42.0,501);  // for 4.4.4.4 blocking 1226 evs
+  //  Chebyshev<CoarseVector>      IRLCheby(0.0025,42.0,501);  // for 4.4.4.4 blocking 1059 evs
+							  //							  3e-4,2);
+  Chebyshev<CoarseVector>      IRLCheby(0.0018,42.0,301);  // for 4.4.4.4 blocking  // 790 evs
+  
   MrhsHermMatrix MrhsCoarseOp     (mrhs);
 
   CoarseVector pm_src(CoarseMrhs);
   pm_src = ComplexD(1.0);
   PowerMethod<CoarseVector>       cPM;   cPM(MrhsCoarseOp,pm_src);
 
-  int Nk=nrhs*30;
+  //  int Nk=nrhs*30; // 4.4.6.4
   //  int Nk=nrhs*80;
-  int Nm=Nk*4;
-  int Nstop=Nk;
+  int Nk=nrhs*60; // 720
+  int Nm=Nk*4;    // 2880 ; generally finishes at 1440
+  int Nstop=512;
   int Nconv_test_interval=1;
   
   ImplicitlyRestartedBlockLanczosCoarse<CoarseVector> IRL(MrhsCoarseOp,
@@ -299,7 +308,7 @@ int main (int argc, char ** argv)
 							  nrhs,
 							  Nk,
 							  Nm,
-							  1e-4,20);
+							  3e-4,2);
 
   std::vector<RealD>            eval(Nm);
   std::vector<CoarseVector>     evec(Nm,Coarse5d);
@@ -331,7 +340,7 @@ int main (int argc, char ** argv)
   // Extra HDCG parameters
   //////////////////////////
   int maxit=3000;
-  ConjugateGradient<CoarseVector>  CG(5.0e-2,maxit,false);
+  ConjugateGradient<CoarseVector>  CG(7.5e-2,maxit,false);
   RealD lo=2.0;
   int ord = 7;