Const correctness patch

Grid/algorithms/LinearOperator.h

@@ -380,12 +380,6 @@ namespace Grid {
     template<class Field> class OperatorFunction {
     public:
       virtual void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) = 0;
-      virtual void operator() (LinearOperatorBase<Field> &Linop, const std::vector<Field> &in,std::vector<Field> &out) {
-	assert(in.size()==out.size());
-	for(int k=0;k<in.size();k++){
-	  (*this)(Linop,in[k],out[k]);
-	}
-      };
     };
 
     template<class Field> class LinearFunction {
@@ -427,7 +421,7 @@ namespace Grid {
   // Hermitian operator Linear function and operator function
   ////////////////////////////////////////////////////////////////////////////////////////////
     template<class Field>
-    class HermOpOperatorFunction : public OperatorFunction<Field> {
+      class HermOpOperatorFunction : public OperatorFunction<Field> {
       void operator() (LinearOperatorBase<Field> &Linop, const Field &in, Field &out) {
 	Linop.HermOp(in,out);
       };

Grid/algorithms/SparseMatrix.h

@@ -55,14 +55,6 @@ namespace Grid {
     template<class Field> class CheckerBoardedSparseMatrixBase : public SparseMatrixBase<Field> {
     public:
       virtual GridBase *RedBlackGrid(void)=0;
-
-      //////////////////////////////////////////////////////////////////////
-      // Query the even even properties to make algorithmic decisions
-      //////////////////////////////////////////////////////////////////////
-      virtual RealD  Mass(void)        { return 0.0; };
-      virtual int    ConstEE(void)     { return 0; }; // Disable assumptions unless overridden
-      virtual int    isTrivialEE(void) { return 0; }; // by a derived class that knows better
-
       // half checkerboard operaions
       virtual  void Meooe    (const Field &in, Field &out)=0;
       virtual  void Mooee    (const Field &in, Field &out)=0;

Grid/algorithms/iterative/BlockConjugateGradient.h

@@ -33,7 +33,7 @@ directory
 
 namespace Grid {
 
-enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec, BlockCGrQVec };
+enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS };
 
 //////////////////////////////////////////////////////////////////////////
 // Block conjugate gradient. Dimension zero should be the block direction
@@ -42,6 +42,7 @@ template <class Field>
 class BlockConjugateGradient : public OperatorFunction<Field> {
  public:
 
+
   typedef typename Field::scalar_type scomplex;
 
   int blockDim ;
@@ -53,15 +54,21 @@ 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),PrintInterval(100)
+    : Tolerance(tol), CGtype(cgtype),   blockDim(_Orthog),  MaxIterations(maxit), ErrorOnNoConverge(err_on_no_conv)
   {};
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Thin QR factorisation (google it)
 ////////////////////////////////////////////////////////////////////////////////////////////////////
+void ThinQRfact (Eigen::MatrixXcd &m_rr,
+		 Eigen::MatrixXcd &C,
+		 Eigen::MatrixXcd &Cinv,
+		 Field & Q,
+		 const Field & R)
+{
+  int Orthog = blockDim; // First dimension is block dim; this is an assumption
   ////////////////////////////////////////////////////////////////////////////////////////////////////
   //Dimensions
   // R_{ferm x Nblock} =  Q_{ferm x Nblock} x  C_{Nblock x Nblock} -> ferm x Nblock
@@ -78,20 +85,22 @@ class BlockConjugateGradient : public OperatorFunction<Field> {
   // Cdag C = Rdag R ; passes.
   // QdagQ  = 1      ; passes
   ////////////////////////////////////////////////////////////////////////////////////////////////////
-void ThinQRfact (Eigen::MatrixXcd &m_rr,
-		 Eigen::MatrixXcd &C,
-		 Eigen::MatrixXcd &Cinv,
-		 Field & Q,
-		 const Field & R)
-{
-  int Orthog = blockDim; // First dimension is block dim; this is an assumption
   sliceInnerProductMatrix(m_rr,R,R,Orthog);
 
   // Force manifest hermitian to avoid rounding related
   m_rr = 0.5*(m_rr+m_rr.adjoint());
 
-  Eigen::MatrixXcd L    = m_rr.llt().matrixL(); 
+#if 0
+  std::cout << " Calling Cholesky  ldlt on m_rr "  << m_rr <<std::endl;
+  Eigen::MatrixXcd L_ldlt = m_rr.ldlt().matrixL(); 
+  std::cout << " Called Cholesky  ldlt on m_rr "  << L_ldlt <<std::endl;
+  auto  D_ldlt = m_rr.ldlt().vectorD(); 
+  std::cout << " Called Cholesky  ldlt on m_rr "  << D_ldlt <<std::endl;
+#endif
 
+  //  std::cout << " Calling Cholesky  llt on m_rr "  <<std::endl;
+  Eigen::MatrixXcd L    = m_rr.llt().matrixL(); 
+  //  std::cout << " Called Cholesky  llt on m_rr "  << L <<std::endl;
   C    = L.adjoint();
   Cinv = C.inverse();
   ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -103,25 +112,6 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr,
   ////////////////////////////////////////////////////////////////////////////////////////////////////
   sliceMulMatrix(Q,Cinv,R,Orthog);
 }
-// see comments above
-void ThinQRfact (Eigen::MatrixXcd &m_rr,
-		 Eigen::MatrixXcd &C,
-		 Eigen::MatrixXcd &Cinv,
-		 std::vector<Field> & Q,
-		 const std::vector<Field> & R)
-{
-  InnerProductMatrix(m_rr,R,R);
-
-  m_rr = 0.5*(m_rr+m_rr.adjoint());
-
-  Eigen::MatrixXcd L    = m_rr.llt().matrixL(); 
-
-  C    = L.adjoint();
-  Cinv = C.inverse();
-
-  MulMatrix(Q,Cinv,R);
-}
-
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // Call one of several implementations
 ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -129,20 +119,14 @@ void operator()(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi)
 {
   if ( CGtype == BlockCGrQ ) {
     BlockCGrQsolve(Linop,Src,Psi);
+  } else if (CGtype == BlockCG ) {
+    BlockCGsolve(Linop,Src,Psi);
   } else if (CGtype == CGmultiRHS ) {
     CGmultiRHSsolve(Linop,Src,Psi);
   } else {
     assert(0);
   }
 }
-virtual void operator()(LinearOperatorBase<Field> &Linop, const std::vector<Field> &Src, std::vector<Field> &Psi) 
-{
-  if ( CGtype == BlockCGrQVec ) {
-    BlockCGrQsolveVec(Linop,Src,Psi);
-  } else {
-    assert(0);
-  }
-}
 
 ////////////////////////////////////////////////////////////////////////////
 // BlockCGrQ implementation:
@@ -155,8 +139,7 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
 {
   int Orthog = blockDim; // First dimension is block dim; this is an assumption
   Nblock = B._grid->_fdimensions[Orthog];
-/* FAKE */
+
-  Nblock=8;
   std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
 
   X.checkerboard = B.checkerboard;
@@ -219,10 +202,15 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   std::cout << GridLogMessage<<"BlockCGrQ algorithm initialisation " <<std::endl;
 
   //1.  QC = R = B-AX, D = Q     ; QC => Thin QR factorisation (google it)
+
   Linop.HermOp(X, AD);
   tmp = B - AD;  
-
+  //std::cout << GridLogMessage << " initial tmp " << norm2(tmp)<< std::endl;
   ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
+  //std::cout << GridLogMessage << " initial Q " << norm2(Q)<< std::endl;
+  //std::cout << GridLogMessage << " m_rr " << m_rr<<std::endl;
+  //std::cout << GridLogMessage << " m_C " << m_C<<std::endl;
+  //std::cout << GridLogMessage << " m_Cinv " << m_Cinv<<std::endl;
   D=Q;
 
   std::cout << GridLogMessage<<"BlockCGrQ computed initial residual and QR fact " <<std::endl;
@@ -244,12 +232,14 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
     MatrixTimer.Start();
     Linop.HermOp(D, Z);      
     MatrixTimer.Stop();
+    //std::cout << GridLogMessage << " norm2 Z " <<norm2(Z)<<std::endl;
 
     //4. M  = [D^dag Z]^{-1}
     sliceInnerTimer.Start();
     sliceInnerProductMatrix(m_DZ,D,Z,Orthog);
     sliceInnerTimer.Stop();
     m_M       = m_DZ.inverse();
+    //std::cout << GridLogMessage << " m_DZ " <<m_DZ<<std::endl;
     
     //5. X  = X + D MC
     m_tmp     = m_M * m_C;
@@ -267,7 +257,6 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
     
     //7. D  = Q + D S^dag
     m_tmp = m_S.adjoint();
-
     sliceMaddTimer.Start();
     sliceMaddMatrix(D,m_tmp,D,Q,Orthog);
     sliceMaddTimer.Stop();
@@ -328,6 +317,152 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   IterationsToComplete = k;
 }
 //////////////////////////////////////////////////////////////////////////
+// Block conjugate gradient; Original O'Leary Dimension zero should be the block direction
+//////////////////////////////////////////////////////////////////////////
+void BlockCGsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &Psi) 
+{
+  int Orthog = blockDim; // First dimension is block dim; this is an assumption
+  Nblock = Src._grid->_fdimensions[Orthog];
+
+  std::cout<<GridLogMessage<<" Block Conjugate Gradient : Orthog "<<Orthog<<" Nblock "<<Nblock<<std::endl;
+
+  Psi.checkerboard = Src.checkerboard;
+  conformable(Psi, Src);
+
+  Field P(Src);
+  Field AP(Src);
+  Field R(Src);
+  
+  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);
+  RealD sssum=0;
+  for(int b=0;b<Nblock;b++) sssum+=ssq[b];
+
+  sliceNorm(residuals,Src,Orthog);
+  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
+
+  sliceNorm(residuals,Psi,Orthog);
+  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
+
+  // Initial search dir is guess
+  Linop.HermOp(Psi, AP);
+  
+
+  /************************************************************************
+   * 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
+   */
+
+  R = Src - AP;  
+  P = R;
+  sliceInnerProductMatrix(m_rr,R,R,Orthog);
+
+  GridStopWatch sliceInnerTimer;
+  GridStopWatch sliceMaddTimer;
+  GridStopWatch MatrixTimer;
+  GridStopWatch SolverTimer;
+  SolverTimer.Start();
+
+  int k;
+  for (k = 1; k <= MaxIterations; k++) {
+
+    RealD rrsum=0;
+    for(int b=0;b<Nblock;b++) rrsum+=real(m_rr(b,b));
+
+    std::cout << GridLogIterative << "\titeration "<<k<<" rr_sum "<<rrsum<<" ssq_sum "<< sssum
+	      <<" / "<<std::sqrt(rrsum/sssum) <<std::endl;
+
+    MatrixTimer.Start();
+    Linop.HermOp(P, AP);
+    MatrixTimer.Stop();
+
+    // Alpha
+    sliceInnerTimer.Start();
+    sliceInnerProductMatrix(m_pAp,P,AP,Orthog);
+    sliceInnerTimer.Stop();
+    m_pAp_inv = m_pAp.inverse();
+    m_alpha   = m_pAp_inv * m_rr ;
+
+    // Psi, R update
+    sliceMaddTimer.Start();
+    sliceMaddMatrix(Psi,m_alpha, P,Psi,Orthog);     // add alpha *  P to psi
+    sliceMaddMatrix(R  ,m_alpha,AP,  R,Orthog,-1.0);// sub alpha * AP to resid
+    sliceMaddTimer.Stop();
+
+    // Beta
+    m_rr_inv = m_rr.inverse();
+    sliceInnerTimer.Start();
+    sliceInnerProductMatrix(m_rr,R,R,Orthog);
+    sliceInnerTimer.Stop();
+    m_beta = m_rr_inv *m_rr;
+
+    // Search update
+    sliceMaddTimer.Start();
+    sliceMaddMatrix(AP,m_beta,P,R,Orthog);
+    sliceMaddTimer.Stop();
+    P= 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;
+
+      Linop.HermOp(Psi, AP);
+      AP = AP-Src;
+      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(norm2(AP)/norm2(Src)) <<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;
+}
+//////////////////////////////////////////////////////////////////////////
 // multiRHS conjugate gradient. Dimension zero should be the block direction
 // Use this for spread out across nodes
 //////////////////////////////////////////////////////////////////////////
@@ -465,233 +600,6 @@ void CGmultiRHSsolve(LinearOperatorBase<Field> &Linop, const Field &Src, Field &
   IterationsToComplete = k;
 }
 
-void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector<Field> &X, const 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]);  
-  }}
-}
-void MaddMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X,const std::vector<Field> &Y,RealD scale=1.0){
-  // Should make this cache friendly with site outermost, parallel_for
-  // Deal with case AP aliases with either Y or X
-  std::vector<Field> tmp(Nblock,X[0]);
-  for(int b=0;b<Nblock;b++){
-    tmp[b]   = Y[b];
-    for(int bp=0;bp<Nblock;bp++) {
-      tmp[b] = tmp[b] + (scale*m(bp,b))*X[bp]; 
-    }
-  }
-  for(int b=0;b<Nblock;b++){
-    AP[b] = tmp[b];
-  }
-}
-void MulMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X){
-  // Should make this cache friendly with site outermost, parallel_for
-  for(int b=0;b<Nblock;b++){
-    AP[b] = zero;
-    for(int bp=0;bp<Nblock;bp++) {
-      AP[b] += (m(bp,b))*X[bp]; 
-    }
-  }
-}
-double normv(const std::vector<Field> &P){
-  double nn = 0.0;
-  for(int b=0;b<Nblock;b++) {
-    nn+=norm2(P[b]);
-  }
-  return nn;
-}
-
-////////////////////////////////////////////////////////////////////////////
-// BlockCGrQvec implementation:
-//--------------------------
-// X is guess/Solution
-// B is RHS
-// Solve A X_i = B_i    ;        i refers to Nblock index
-////////////////////////////////////////////////////////////////////////////
-void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field> &B, std::vector<Field> &X) 
-{
-  Nblock = B.size();
-  assert(Nblock == X.size());
-
-  std::cout<<GridLogMessage<<" Block Conjugate Gradient Vec rQ : Nblock "<<Nblock<<std::endl;
-
-  for(int b=0;b<Nblock;b++){ 
-    X[b].checkerboard = B[b].checkerboard;
-    conformable(X[b], B[b]);
-    conformable(X[b], X[0]); 
-  }
-
-  Field Fake(B[0]);
-
-  std::vector<Field> tmp(Nblock,Fake);
-  std::vector<Field>   Q(Nblock,Fake);
-  std::vector<Field>   D(Nblock,Fake);
-  std::vector<Field>   Z(Nblock,Fake);
-  std::vector<Field>  AD(Nblock,Fake);
-
-  Eigen::MatrixXcd m_DZ     = Eigen::MatrixXcd::Identity(Nblock,Nblock);
-  Eigen::MatrixXcd m_M      = Eigen::MatrixXcd::Identity(Nblock,Nblock);
-  Eigen::MatrixXcd m_rr     = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-
-  Eigen::MatrixXcd m_C      = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-  Eigen::MatrixXcd m_Cinv   = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-  Eigen::MatrixXcd m_S      = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-  Eigen::MatrixXcd m_Sinv   = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-
-  Eigen::MatrixXcd m_tmp    = Eigen::MatrixXcd::Identity(Nblock,Nblock);
-  Eigen::MatrixXcd m_tmp1   = Eigen::MatrixXcd::Identity(Nblock,Nblock);
-
-  // Initial residual computation & set up
-  std::vector<RealD> residuals(Nblock);
-  std::vector<RealD> ssq(Nblock);
-
-  RealD sssum=0;
-  for(int b=0;b<Nblock;b++){ ssq[b] = norm2(B[b]);}
-  for(int b=0;b<Nblock;b++) sssum+=ssq[b];
-
-  for(int b=0;b<Nblock;b++){ residuals[b] = norm2(B[b]);}
-  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
-
-  for(int b=0;b<Nblock;b++){ residuals[b] = norm2(X[b]);}
-  for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
-
-  /************************************************************************
-   * Block conjugate gradient rQ (Sebastien Birk Thesis, after Dubrulle 2001)
-   ************************************************************************
-   * Dimensions:
-   *
-   *   X,B==(Nferm x Nblock)
-   *   A==(Nferm x Nferm)
-   *  
-   * Nferm = Nspin x Ncolour x Ncomplex x Nlattice_site
-   * 
-   * QC = R = B-AX, D = Q     ; QC => Thin QR factorisation (google it)
-   * for k: 
-   *   Z  = AD
-   *   M  = [D^dag Z]^{-1}
-   *   X  = X + D MC
-   *   QS = Q - ZM
-   *   D  = Q + D S^dag
-   *   C  = S C
-   */
-  ///////////////////////////////////////
-  // Initial block: initial search dir is guess
-  ///////////////////////////////////////
-  std::cout << GridLogMessage<<"BlockCGrQvec algorithm initialisation " <<std::endl;
-
-  //1.  QC = R = B-AX, D = Q     ; QC => Thin QR factorisation (google it)
-  for(int b=0;b<Nblock;b++) {
-    Linop.HermOp(X[b], AD[b]);
-    tmp[b] = B[b] - AD[b];  
-  }
-
-  ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
-
-  for(int b=0;b<Nblock;b++) D[b]=Q[b];
-
-  std::cout << GridLogMessage<<"BlockCGrQ vec computed initial residual and QR fact " <<std::endl;
-
-  ///////////////////////////////////////
-  // Timers
-  ///////////////////////////////////////
-  GridStopWatch sliceInnerTimer;
-  GridStopWatch sliceMaddTimer;
-  GridStopWatch QRTimer;
-  GridStopWatch MatrixTimer;
-  GridStopWatch SolverTimer;
-  SolverTimer.Start();
-
-  int k;
-  for (k = 1; k <= MaxIterations; k++) {
-
-    //3. Z  = AD
-    MatrixTimer.Start();
-    for(int b=0;b<Nblock;b++) Linop.HermOp(D[b], Z[b]);      
-    MatrixTimer.Stop();
-
-    //4. M  = [D^dag Z]^{-1}
-    sliceInnerTimer.Start();
-    InnerProductMatrix(m_DZ,D,Z);
-    sliceInnerTimer.Stop();
-    m_M       = m_DZ.inverse();
-    
-    //5. X  = X + D MC
-    m_tmp     = m_M * m_C;
-    sliceMaddTimer.Start();
-    MaddMatrix(X,m_tmp, D,X);     
-    sliceMaddTimer.Stop();
-
-    //6. QS = Q - ZM
-    sliceMaddTimer.Start();
-    MaddMatrix(tmp,m_M,Z,Q,-1.0);
-    sliceMaddTimer.Stop();
-    QRTimer.Start();
-    ThinQRfact (m_rr, m_S, m_Sinv, Q, tmp);
-    QRTimer.Stop();
-    
-    //7. D  = Q + D S^dag
-    m_tmp = m_S.adjoint();
-    sliceMaddTimer.Start();
-    MaddMatrix(D,m_tmp,D,Q);
-    sliceMaddTimer.Stop();
-
-    //8. C  = S C
-    m_C = m_S*m_C;
-    
-    /*********************
-     * convergence monitor
-     *********************
-     */
-    m_rr = m_C.adjoint() * m_C;
-
-    RealD max_resid=0;
-    RealD rrsum=0;
-    RealD rr;
-
-    for(int b=0;b<Nblock;b++) {
-      rrsum+=real(m_rr(b,b));
-      rr = real(m_rr(b,b))/ssq[b];
-      if ( rr > max_resid ) max_resid = rr;
-    }
-
-    std::cout << GridLogIterative << "\t Block Iteration "<<k<<" ave resid "<< sqrt(rrsum/sssum) << " max "<< sqrt(max_resid) <<std::endl;
-
-    if ( max_resid < Tolerance*Tolerance ) { 
-
-      SolverTimer.Stop();
-
-      std::cout << GridLogMessage<<"BlockCGrQ 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(X[b], AD[b]);
-      for(int b=0;b<Nblock;b++) AD[b] = AD[b]-B[b];
-      std::cout << GridLogMessage <<"\t True residual is " << std::sqrt(normv(AD)/normv(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;
-      std::cout << GridLogMessage << "\tThinQRfact " << QRTimer.Elapsed()  <<std::endl;
-	    
-      IterationsToComplete = k;
-      return;
-    }
-
-  }
-  std::cout << GridLogMessage << "BlockConjugateGradient(rQ) did NOT converge" << std::endl;
-
-  if (ErrorOnNoConverge) assert(0);
-  IterationsToComplete = k;
-}
-
-
-
 };
 
 }

Grid/algorithms/iterative/ConjugateGradient.h

@@ -133,7 +133,7 @@ class ConjugateGradient : public OperatorFunction<Field> {
       LinalgTimer.Stop();
 
       std::cout << GridLogIterative << "ConjugateGradient: Iteration " << k
-                << " residual^2 " << sqrt(cp/ssq) << " target " << Tolerance << std::endl;
+                << " residual " << cp << " target " << rsq << std::endl;
 
       // Stopping condition
       if (cp <= rsq) {
@@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
 	std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
 	std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
 
-        std::cout << GridLogMessage << "Time breakdown "<<std::endl;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
-	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 

Grid/algorithms/iterative/SchurRedBlack.h

@@ -86,23 +86,229 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
    */
 namespace Grid {
 
-  ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  // Use base class to share code
-  ///////////////////////////////////////////////////////////////////////////////////////////////////////
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
   // Take a matrix and form a Red Black solver calling a Herm solver
   // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  template<class Field> class SchurRedBlackBase {
+  // Now make the norm reflect extra factor of Mee
-  protected:
+  template<class Field> class SchurRedBlackStaggeredSolve {
-    typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
+  private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
     bool subGuess;
   public:
 
-    SchurRedBlackBase(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
+    /////////////////////////////////////////////////////
-    _HermitianRBSolver(HermitianRBSolver) 
+    // Wrap the usual normal equations Schur trick
+    /////////////////////////////////////////////////////
+  SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
+     _HermitianRBSolver(HermitianRBSolver) 
+    { 
+      CBfactorise=0;
+      subtractGuess(initSubGuess);
+    };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
+
+    template<class Matrix>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out){
+      ZeroGuesser<Field> guess;
+      (*this)(_Matrix,in,out,guess);
+    }
+    template<class Matrix, class Guesser>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out, Guesser &guess){
+
+      // FIXME CGdiagonalMee not implemented virtual function
+      // FIXME use CBfactorise to control schur decomp
+      GridBase *grid = _Matrix.RedBlackGrid();
+      GridBase *fgrid= _Matrix.Grid();
+
+      SchurStaggeredOperator<Matrix,Field> _HermOpEO(_Matrix);
+ 
+      Field src_e(grid);
+      Field src_o(grid);
+      Field sol_e(grid);
+      Field sol_o(grid);
+      Field   tmp(grid);
+      Field  Mtmp(grid);
+      Field resid(fgrid);
+      
+      std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve " <<std::endl;
+      pickCheckerboard(Even,src_e,in);
+      pickCheckerboard(Odd ,src_o,in);
+      pickCheckerboard(Even,sol_e,out);
+      pickCheckerboard(Odd ,sol_o,out);
+      std::cout << GridLogMessage << " SchurRedBlackStaggeredSolve checkerboards picked" <<std::endl;
+    
+      /////////////////////////////////////////////////////
+      // src_o = (source_o - Moe MeeInv source_e)
+      /////////////////////////////////////////////////////
+      _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.checkerboard ==Even);
+      _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
+      tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
+
+      //src_o = tmp;     assert(src_o.checkerboard ==Odd);
+      _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm.
+
+      //////////////////////////////////////////////////////////////
+      // Call the red-black solver
+      //////////////////////////////////////////////////////////////
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
+      guess(src_o, sol_o);
+      Mtmp = sol_o;
+      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called  the Mpc solver" <<std::endl;
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
+
+      ///////////////////////////////////////////////////
+      // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+      ///////////////////////////////////////////////////
+      _Matrix.Meooe(sol_o,tmp);        assert(  tmp.checkerboard   ==Even);
+      src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
+      _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
+     
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver reconstructed other CB" <<std::endl;
+      setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
+      setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
+      std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl;
+
+      // Verify the unprec residual
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
+        std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
+    }     
+  };
+  template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Take a matrix and form a Red Black solver calling a Herm solver
+  // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  template<class Field> class SchurRedBlackDiagMooeeSolve {
+  private:
+    OperatorFunction<Field> & _HermitianRBSolver;
+    int CBfactorise;
+    bool subGuess;
+  public:
+
+    /////////////////////////////////////////////////////
+    // Wrap the usual normal equations Schur trick
+    /////////////////////////////////////////////////////
+  SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0, const bool initSubGuess = false)  :  _HermitianRBSolver(HermitianRBSolver) 
+  { 
+    CBfactorise=cb;
+    subtractGuess(initSubGuess);
+  };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
+    template<class Matrix>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out){
+      ZeroGuesser<Field> guess;
+      (*this)(_Matrix,in,out,guess);
+    }
+    template<class Matrix, class Guesser>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){
+
+      // FIXME CGdiagonalMee not implemented virtual function
+      // FIXME use CBfactorise to control schur decomp
+      GridBase *grid = _Matrix.RedBlackGrid();
+      GridBase *fgrid= _Matrix.Grid();
+
+      SchurDiagMooeeOperator<Matrix,Field> _HermOpEO(_Matrix);
+ 
+      Field src_e(grid);
+      Field src_o(grid);
+      Field sol_e(grid);
+      Field sol_o(grid);
+      Field   tmp(grid);
+      Field  Mtmp(grid);
+      Field resid(fgrid);
+
+      pickCheckerboard(Even,src_e,in);
+      pickCheckerboard(Odd ,src_o,in);
+      pickCheckerboard(Even,sol_e,out);
+      pickCheckerboard(Odd ,sol_o,out);
+    
+      /////////////////////////////////////////////////////
+      // src_o = Mdag * (source_o - Moe MeeInv source_e)
+      /////////////////////////////////////////////////////
+      _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.checkerboard ==Even);
+      _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
+      tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
+
+      // get the right MpcDag
+      _HermOpEO.MpcDag(tmp,src_o);     assert(src_o.checkerboard ==Odd);       
+
+      //////////////////////////////////////////////////////////////
+      // Call the red-black solver
+      //////////////////////////////////////////////////////////////
+      std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
+      guess(src_o,sol_o);
+      Mtmp = sol_o;
+      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
+
+      ///////////////////////////////////////////////////
+      // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
+      ///////////////////////////////////////////////////
+      _Matrix.Meooe(sol_o,tmp);        assert(  tmp.checkerboard   ==Even);
+      src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
+      _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
+     
+      setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
+      setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
+
+      // Verify the unprec residual
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
+
+        std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
+    }     
+  };
+
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Take a matrix and form a Red Black solver calling a Herm solver
+  // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////
+  template<class Field> class SchurRedBlackDiagTwoSolve {
+  private:
+    OperatorFunction<Field> & _HermitianRBSolver;
+    int CBfactorise;
+    bool subGuess;
+  public:
+
+    /////////////////////////////////////////////////////
+    // Wrap the usual normal equations Schur trick
+    /////////////////////////////////////////////////////
+  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
+     _HermitianRBSolver(HermitianRBSolver) 
     { 
       CBfactorise = 0;
       subtractGuess(initSubGuess);
@@ -116,86 +322,12 @@ namespace Grid {
       return subGuess;
     }
 
-    /////////////////////////////////////////////////////////////
+    template<class Matrix>
-    // Shared code
-    /////////////////////////////////////////////////////////////
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
       ZeroGuesser<Field> guess;
       (*this)(_Matrix,in,out,guess);
     }
-    void operator()(Matrix &_Matrix, const std::vector<Field> &in, std::vector<Field> &out) 
+    template<class Matrix,class Guesser>
-    {
-      ZeroGuesser<Field> guess;
-      (*this)(_Matrix,in,out,guess);
-    }
-
-    template<class Guesser>
-    void operator()(Matrix &_Matrix, const std::vector<Field> &in, std::vector<Field> &out,Guesser &guess) 
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      GridBase *fgrid= _Matrix.Grid();
-      int nblock = in.size();
-
-      std::vector<Field> src_o(nblock,grid);
-      std::vector<Field> sol_o(nblock,grid);
-      
-      std::vector<Field> guess_save;
-
-      Field resid(fgrid);
-      Field tmp(grid);
-
-      ////////////////////////////////////////////////
-      // Prepare RedBlack source
-      ////////////////////////////////////////////////
-      for(int b=0;b<nblock;b++){
-	RedBlackSource(_Matrix,in[b],tmp,src_o[b]);
-      }
-      ////////////////////////////////////////////////
-      // Make the guesses
-      ////////////////////////////////////////////////
-      if ( subGuess ) guess_save.resize(nblock,grid);
-
-      for(int b=0;b<nblock;b++){
-	guess(src_o[b],sol_o[b]); 
-
-	if ( subGuess ) { 
-	  guess_save[b] = sol_o[b];
-	}
-      }
-      //////////////////////////////////////////////////////////////
-      // Call the block solver
-      //////////////////////////////////////////////////////////////
-      std::cout<<GridLogMessage << "SchurRedBlackBase calling the solver for "<<nblock<<" RHS" <<std::endl;
-      RedBlackSolve(_Matrix,src_o,sol_o);
-
-      ////////////////////////////////////////////////
-      // A2A boolean behavioural control & reconstruct other checkerboard
-      ////////////////////////////////////////////////
-      for(int b=0;b<nblock;b++) {
-
-	if (subGuess)   sol_o[b] = sol_o[b] - guess_save[b];
-
-	///////// Needs even source //////////////
-	pickCheckerboard(Even,tmp,in[b]);
-	RedBlackSolution(_Matrix,sol_o[b],tmp,out[b]);
-
-	/////////////////////////////////////////////////
-	// Check unprec residual if possible
-	/////////////////////////////////////////////////
-	if ( ! subGuess ) {
-	  _Matrix.M(out[b],resid); 
-	  resid = resid-in[b];
-	  RealD ns = norm2(in[b]);
-	  RealD nr = norm2(resid);
-	
-	  std::cout<<GridLogMessage<< "SchurRedBlackBase solver true unprec resid["<<b<<"] "<<std::sqrt(nr/ns) << std::endl;
-	} else {
-	  std::cout<<GridLogMessage<< "SchurRedBlackBase Guess subtracted after solve["<<b<<"] " << std::endl;
-	}
-
-      }
-    }
-    template<class Guesser>
     void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){
 
       // FIXME CGdiagonalMee not implemented virtual function
@@ -203,105 +335,42 @@ namespace Grid {
       GridBase *grid = _Matrix.RedBlackGrid();
       GridBase *fgrid= _Matrix.Grid();
 
-      Field resid(fgrid);
+      SchurDiagTwoOperator<Matrix,Field> _HermOpEO(_Matrix);
-      Field src_o(grid);
+ 
       Field src_e(grid);
+      Field src_o(grid);
+      Field sol_e(grid);
       Field sol_o(grid);
-
-      ////////////////////////////////////////////////
-      // RedBlack source
-      ////////////////////////////////////////////////
-      RedBlackSource(_Matrix,in,src_e,src_o);
-
-      ////////////////////////////////
-      // Construct the guess
-      ////////////////////////////////
-      Field   tmp(grid);
-      guess(src_o,sol_o);
-
-      Field  guess_save(grid);
-      guess_save = sol_o;
-
-      //////////////////////////////////////////////////////////////
-      // Call the red-black solver
-      //////////////////////////////////////////////////////////////
-      RedBlackSolve(_Matrix,src_o,sol_o);
-
-      ////////////////////////////////////////////////
-      // Fionn A2A boolean behavioural control
-      ////////////////////////////////////////////////
-      if (subGuess)      sol_o= sol_o-guess_save;
-
-      ///////////////////////////////////////////////////
-      // RedBlack solution needs the even source
-      ///////////////////////////////////////////////////
-      RedBlackSolution(_Matrix,sol_o,src_e,out);
-
-      // Verify the unprec residual
-      if ( ! subGuess ) {
-        _Matrix.M(out,resid); 
-        resid = resid-in;
-        RealD ns = norm2(in);
-        RealD nr = norm2(resid);
-
-        std::cout<<GridLogMessage << "SchurRedBlackBase solver true unprec resid "<< std::sqrt(nr/ns) << std::endl;
-      } else {
-        std::cout << GridLogMessage << "SchurRedBlackBase Guess subtracted after solve." << std::endl;
-      }
-    }     
-    
-    /////////////////////////////////////////////////////////////
-    // Override in derived. Not virtual as template methods
-    /////////////////////////////////////////////////////////////
-    virtual void RedBlackSource  (Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)                =0;
-    virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol)          =0;
-    virtual void RedBlackSolve   (Matrix & _Matrix,const Field &src_o, Field &sol_o)                           =0;
-    virtual void RedBlackSolve   (Matrix & _Matrix,const std::vector<Field> &src_o,  std::vector<Field> &sol_o)=0;
-
-  };
-
-  template<class Field> class SchurRedBlackStaggeredSolve : public SchurRedBlackBase<Field> {
-  public:
-    typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
-
-    SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false) 
-      :    SchurRedBlackBase<Field> (HermitianRBSolver,initSubGuess) 
-    {
-    }
-
-    //////////////////////////////////////////////////////
-    // Override RedBlack specialisation
-    //////////////////////////////////////////////////////
-    virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      GridBase *fgrid= _Matrix.Grid();
-
       Field   tmp(grid);
       Field  Mtmp(grid);
+      Field resid(fgrid);
 
-      pickCheckerboard(Even,src_e,src);
+      pickCheckerboard(Even,src_e,in);
-      pickCheckerboard(Odd ,src_o,src);
+      pickCheckerboard(Odd ,src_o,in);
+      pickCheckerboard(Even,sol_e,out);
+      pickCheckerboard(Odd ,sol_o,out);
     
       /////////////////////////////////////////////////////
-      // src_o = (source_o - Moe MeeInv source_e)
+      // src_o = Mdag * (source_o - Moe MeeInv source_e)
       /////////////////////////////////////////////////////
       _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.checkerboard ==Even);
       _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
       tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
 
-      _Matrix.Mooee(tmp,src_o); // Extra factor of "m" in source from dumb choice of matrix norm.
+      // 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_c,Field &sol)
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      GridBase *fgrid= _Matrix.Grid();
 
-      Field   tmp(grid);
+      //////////////////////////////////////////////////////////////
-      Field   sol_e(grid);
+      // Call the red-black solver
-      Field   src_e(grid);
+      //////////////////////////////////////////////////////////////
-
+      std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
-      src_e = src_e_c; // Const correctness
+//      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      guess(src_o,tmp);
+      Mtmp = tmp;
+      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)      tmp = tmp-Mtmp;
+      _Matrix.MooeeInv(tmp,sol_o);       assert(  sol_o.checkerboard   ==Odd);
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -310,116 +379,78 @@ namespace Grid {
       src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
       _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
      
-      setCheckerboard(sol,sol_e); assert(  sol_e.checkerboard ==Even);
+      setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
-      setCheckerboard(sol,sol_o); assert(  sol_o.checkerboard ==Odd );
+      setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
-    }
+
-    virtual void RedBlackSolve   (Matrix & _Matrix,const Field &src_o, Field &sol_o)
+      // Verify the unprec residual
-    {
+      if ( ! subGuess ) {
-      SchurStaggeredOperator<Matrix,Field> _HermOpEO(_Matrix);
+        _Matrix.M(out,resid); 
-      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+        resid = resid-in;
-    };
+        RealD ns = norm2(in);
-    virtual void RedBlackSolve   (Matrix & _Matrix,const std::vector<Field> &src_o,  std::vector<Field> &sol_o)
+        RealD nr = norm2(resid);
-    {
+
-      SchurStaggeredOperator<Matrix,Field> _HermOpEO(_Matrix);
+        std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
-      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); 
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
-  template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
-
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  // Site diagonal has Mooee on it.
+  // Take a matrix and form a Red Black solver calling a Herm solver
+  // Use of RB info prevents making SchurRedBlackSolve conform to standard interface
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  template<class Field> class SchurRedBlackDiagMooeeSolve : public SchurRedBlackBase<Field> {
+  template<class Field> class SchurRedBlackDiagTwoMixed {
+  private:
+    LinearFunction<Field> & _HermitianRBSolver;
+    int CBfactorise;
+    bool subGuess;
   public:
-    typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
-
-    SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  
-      : SchurRedBlackBase<Field> (HermitianRBSolver,initSubGuess) {};
-
-
-    //////////////////////////////////////////////////////
-    // Override RedBlack specialisation
-    //////////////////////////////////////////////////////
-    virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)
-    {
-      GridBase *grid = _Matrix.RedBlackGrid();
-      GridBase *fgrid= _Matrix.Grid();
-
-      Field   tmp(grid);
-      Field  Mtmp(grid);
-
-      pickCheckerboard(Even,src_e,src);
-      pickCheckerboard(Odd ,src_o,src);
-
-      /////////////////////////////////////////////////////
-      // src_o = Mdag * (source_o - Moe MeeInv source_e)
-      /////////////////////////////////////////////////////
-      _Matrix.MooeeInv(src_e,tmp);     assert(  tmp.checkerboard ==Even);
-      _Matrix.Meooe   (tmp,Mtmp);      assert( Mtmp.checkerboard ==Odd);     
-      tmp=src_o-Mtmp;                  assert(  tmp.checkerboard ==Odd);     
-
-      // get the right MpcDag
-      SchurDiagMooeeOperator<Matrix,Field> _HermOpEO(_Matrix);
-      _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);
-      Field  src_e_i(grid);
-      ///////////////////////////////////////////////////
-      // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
-      ///////////////////////////////////////////////////
-      _Matrix.Meooe(sol_o,tmp);          assert(  tmp.checkerboard   ==Even);
-      src_e_i = src_e-tmp;               assert(  src_e_i.checkerboard ==Even);
-      _Matrix.MooeeInv(src_e_i,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)
-    {
-      SchurDiagMooeeOperator<Matrix,Field> _HermOpEO(_Matrix);
-      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
-    };
-    virtual void RedBlackSolve   (Matrix & _Matrix,const std::vector<Field> &src_o,  std::vector<Field> &sol_o)
-    {
-      SchurDiagMooeeOperator<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
-  //=> psi = MeeInv phi
-  ///////////////////////////////////////////////////////////////////////////////////////////////////////
-  template<class Field> class SchurRedBlackDiagTwoSolve : public SchurRedBlackBase<Field> {
-  public:
-    typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  
+  SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
-    : SchurRedBlackBase<Field>(HermitianRBSolver,initSubGuess) {};
+     _HermitianRBSolver(HermitianRBSolver) 
-
-    virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)
     { 
+      CBfactorise=0;
+      subtractGuess(initSubGuess);
+    };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
+
+    template<class Matrix>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out){
+      ZeroGuesser<Field> guess;
+      (*this)(_Matrix,in,out,guess);
+    }
+    template<class Matrix, class Guesser>
+    void operator() (Matrix & _Matrix,const Field &in, Field &out,Guesser &guess){
+
+      // FIXME CGdiagonalMee not implemented virtual function
+      // FIXME use CBfactorise to control schur decomp
       GridBase *grid = _Matrix.RedBlackGrid();
       GridBase *fgrid= _Matrix.Grid();
 
       SchurDiagTwoOperator<Matrix,Field> _HermOpEO(_Matrix);
  
+      Field src_e(grid);
+      Field src_o(grid);
+      Field sol_e(grid);
+      Field sol_o(grid);
       Field   tmp(grid);
       Field  Mtmp(grid);
+      Field resid(fgrid);
 
-      pickCheckerboard(Even,src_e,src);
+      pickCheckerboard(Even,src_e,in);
-      pickCheckerboard(Odd ,src_o,src);
+      pickCheckerboard(Odd ,src_o,in);
+      pickCheckerboard(Even,sol_e,out);
+      pickCheckerboard(Odd ,sol_o,out);
     
       /////////////////////////////////////////////////////
       // src_o = Mdag * (source_o - Moe MeeInv source_e)
@@ -430,44 +461,43 @@ namespace Grid {
 
       // 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)
+      //////////////////////////////////////////////////////////////
-    {
+      // Call the red-black solver
-      GridBase *grid = _Matrix.RedBlackGrid();
+      //////////////////////////////////////////////////////////////
-      GridBase *fgrid= _Matrix.Grid();
+      std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
-
+//      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
-      Field   sol_o_i(grid);
+//      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
-      Field   tmp(grid);
+      guess(src_o,tmp);
-      Field   sol_e(grid);
+      Mtmp = tmp;
-
+      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
-      ////////////////////////////////////////////////
+      // Fionn A2A boolean behavioural control
-      // MooeeInv due to pecond
+      if (subGuess)      tmp = tmp-Mtmp;
-      ////////////////////////////////////////////////
+      _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
-      _Matrix.MooeeInv(sol_o,tmp);
-      sol_o_i = tmp;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
       ///////////////////////////////////////////////////
-      _Matrix.Meooe(sol_o_i,tmp);    assert(  tmp.checkerboard   ==Even);
+      _Matrix.Meooe(sol_o,tmp);        assert(  tmp.checkerboard   ==Even);
-      tmp = src_e-tmp;               assert(  src_e.checkerboard ==Even);
+      src_e = src_e-tmp;               assert(  src_e.checkerboard ==Even);
-      _Matrix.MooeeInv(tmp,sol_e);   assert(  sol_e.checkerboard ==Even);
+      _Matrix.MooeeInv(src_e,sol_e);   assert(  sol_e.checkerboard ==Even);
      
-      setCheckerboard(sol,sol_e);    assert(  sol_e.checkerboard ==Even);
+      setCheckerboard(out,sol_e); assert(  sol_e.checkerboard ==Even);
-      setCheckerboard(sol,sol_o_i);  assert(  sol_o_i.checkerboard ==Odd );
+      setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
-    };
 
-    virtual void RedBlackSolve   (Matrix & _Matrix,const Field &src_o, Field &sol_o)
+      // Verify the unprec residual
-    {
+      if ( ! subGuess ) {
-      SchurDiagTwoOperator<Matrix,Field> _HermOpEO(_Matrix);
+        _Matrix.M(out,resid); 
-      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);
+        resid = resid-in;
-    };
+        RealD ns = norm2(in);
-    virtual void RedBlackSolve   (Matrix & _Matrix,const std::vector<Field> &src_o,  std::vector<Field> &sol_o)
+        RealD nr = norm2(resid);
-    {
+
-      SchurDiagTwoOperator<Matrix,Field> _HermOpEO(_Matrix);
+        std::cout << GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid " << std::sqrt(nr / ns) << " nr " << nr << " ns " << ns << std::endl;
-      this->_HermitianRBSolver(_HermOpEO,src_o,sol_o); 
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
+
 }
 #endif

Grid/communicator/Communicator_mpi3.cc

@@ -50,6 +50,8 @@ void CartesianCommunicator::Init(int *argc, char ***argv)
       assert(0);
   }
 
+  Grid_quiesce_nodes();
+
   // Never clean up as done once.
   MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world);
 
@@ -122,8 +124,10 @@ CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors,
   // split the communicator
   //////////////////////////////////////////////////////////////////////////////////////////////////////
   //  int Nparent = parent._processors ; 
+  //  std::cout << " splitting from communicator "<<parent.communicator <<std::endl;
   int Nparent;
   MPI_Comm_size(parent.communicator,&Nparent);
+  //  std::cout << " Parent size  "<<Nparent <<std::endl;
 
   int childsize=1;
   for(int d=0;d<processors.size();d++) {
@@ -132,6 +136,8 @@ CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors,
   int Nchild = Nparent/childsize;
   assert (childsize * Nchild == Nparent);
 
+  //  std::cout << " child size  "<<childsize <<std::endl;
+
   std::vector<int> ccoor(_ndimension); // coor within subcommunicator
   std::vector<int> scoor(_ndimension); // coor of split within parent
   std::vector<int> ssize(_ndimension); // coor of split within parent

Grid/communicator/SharedMemoryMPI.cc

@@ -413,7 +413,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     assert(((uint64_t)ptr&0x3F)==0);
     close(fd);
     WorldShmCommBufs[r] =ptr;
-    //    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
   }
   _ShmAlloc=1;
   _ShmAllocBytes  = bytes;
@@ -455,7 +455,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     assert(((uint64_t)ptr&0x3F)==0);
     close(fd);
     WorldShmCommBufs[r] =ptr;
-    //    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+    std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
   }
   _ShmAlloc=1;
   _ShmAllocBytes  = bytes;
@@ -499,7 +499,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 #endif
       void * ptr =  mmap(NULL,size, PROT_READ | PROT_WRITE, mmap_flag, fd, 0);
       
-      //      std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< size<< "bytes)"<<std::endl;
+      std::cout << "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< size<< "bytes)"<<std::endl;
       if ( ptr == (void * )MAP_FAILED ) {       
 	perror("failed mmap");     
 	assert(0);    

Grid/lattice/Lattice_transfer.h

@@ -464,10 +464,8 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
   assert(orthog>=0);
 
   for(int d=0;d<nh;d++){
-    if ( d!=orthog ) {
+    assert(lg->_processors[d]  == hg->_processors[d]);
-      assert(lg->_processors[d]  == hg->_processors[d]);
+    assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-      assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-    }
   }
 
   // the above should guarantee that the operations are local
@@ -487,7 +485,7 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
 
 
 template<class vobj>
-void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
+void ExtractSliceLocal(Lattice<vobj> &lowDim, const Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
 {
   typedef typename vobj::scalar_object sobj;
 
@@ -501,10 +499,8 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slic
   assert(orthog>=0);
 
   for(int d=0;d<nh;d++){
-    if ( d!=orthog ) {
+    assert(lg->_processors[d]  == hg->_processors[d]);
-      assert(lg->_processors[d]  == hg->_processors[d]);
+    assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-      assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-    }
   }
 
   // the above should guarantee that the operations are local
@@ -524,7 +520,7 @@ void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slic
 
 
 template<class vobj>
-void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
+void Replicate(const Lattice<vobj> &coarse,Lattice<vobj> & fine)
 {
   typedef typename vobj::scalar_object sobj;
 

Grid/log/Log.h

@@ -146,11 +146,9 @@ public:
       if ( log.timestamp ) {
 	log.StopWatch->Stop();
 	GridTime now = log.StopWatch->Elapsed();
-	
 	if ( log.timing_mode==1 ) log.StopWatch->Reset();
 	log.StopWatch->Start();
-	stream << log.evidence()
+	stream << log.evidence()<< std::setw(6)<<now << log.background() << " : " ;
-	       << now	       << log.background() << " : " ;
       }
       stream << log.colour();
       return stream;

Grid/parallelIO/IldgIO.h

@@ -233,8 +233,7 @@ class GridLimeReader : public BinaryIO {
 	//	std::cout << " ReadLatticeObject from offset "<<offset << std::endl;
 	BinarySimpleMunger<sobj,sobj> munge;
 	BinaryIO::readLatticeObject< vobj, sobj >(field, filename, munge, offset, format,nersc_csum,scidac_csuma,scidac_csumb);
-  std::cout << GridLogMessage << "SciDAC checksum A " << std::hex << scidac_csuma << std::dec << std::endl;
+
-  std::cout << GridLogMessage << "SciDAC checksum B " << std::hex << scidac_csumb << std::dec << std::endl;
 	/////////////////////////////////////////////
 	// Insist checksum is next record
 	/////////////////////////////////////////////

Grid/perfmon/Timer.h

@@ -49,34 +49,20 @@ inline double usecond(void) {
 
 typedef  std::chrono::system_clock          GridClock;
 typedef  std::chrono::time_point<GridClock> GridTimePoint;
-
-typedef  std::chrono::seconds               GridSecs;
 typedef  std::chrono::milliseconds          GridMillisecs;
-typedef  std::chrono::microseconds          GridUsecs;
 typedef  std::chrono::microseconds          GridTime;
+typedef  std::chrono::microseconds          GridUsecs;
 
-inline std::ostream& operator<< (std::ostream & stream, const GridSecs & time)
+inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milliseconds & time)
 {
-  stream << time.count()<<" s";
+  stream << time.count()<<" ms";
   return stream;
 }
-inline std::ostream& operator<< (std::ostream & stream, const GridMillisecs & now)
+inline std::ostream& operator<< (std::ostream & stream, const std::chrono::microseconds & time)
 {
-  GridSecs second(1);
+  stream << time.count()<<" usec";
-  auto     secs       = now/second ; 
-  auto     subseconds = now%second ; 
-  stream << secs<<"."<<std::setw(3)<<std::setfill('0')<<subseconds.count()<<" s";
   return stream;
 }
-inline std::ostream& operator<< (std::ostream & stream, const GridUsecs & now)
-{
-  GridSecs second(1);
-  auto     seconds    = now/second ; 
-  auto     subseconds = now%second ; 
-  stream << seconds<<"."<<std::setw(6)<<std::setfill('0')<<subseconds.count()<<" s";
-  return stream;
-}
-
  
 class GridStopWatch {
 private:

Grid/qcd/action/ActionParams.h

@@ -44,15 +44,12 @@ namespace QCD {
   
   struct WilsonImplParams {
     bool overlapCommsCompute;
-    std::vector<Real> twist_n_2pi_L;
     std::vector<Complex> boundary_phases;
     WilsonImplParams() : overlapCommsCompute(false) {
       boundary_phases.resize(Nd, 1.0);
-      twist_n_2pi_L.resize(Nd, 0.0);
     };
-    WilsonImplParams(const std::vector<Complex> phi) : boundary_phases(phi), overlapCommsCompute(false) {
+    WilsonImplParams(const std::vector<Complex> phi)
-      twist_n_2pi_L.resize(Nd, 0.0);
+      : boundary_phases(phi), overlapCommsCompute(false) {}
-    }
   };
 
   struct StaggeredImplParams {

Grid/qcd/action/fermion/FermionOperator.h

@@ -64,6 +64,11 @@ namespace Grid {
       virtual RealD  M    (const FermionField &in, FermionField &out)=0;
       virtual RealD  Mdag (const FermionField &in, FermionField &out)=0;
 
+      // Query the even even properties to make algorithmic decisions
+      virtual int    ConstEE(void) { return 1; }; // clover returns zero as EE depends on gauge field
+      virtual int    isTrivialEE(void) { return 0; };
+      virtual RealD  Mass(void) {return 0.0;};
+
       // half checkerboard operaions
       virtual void   Meooe       (const FermionField &in, FermionField &out)=0;
       virtual void   MeooeDag    (const FermionField &in, FermionField &out)=0;

Grid/qcd/action/fermion/FermionOperatorImpl.h

@@ -240,30 +240,16 @@ namespace QCD {
       GaugeLinkField tmp(GaugeGrid);
 
       Lattice<iScalar<vInteger> > coor(GaugeGrid);
-      ////////////////////////////////////////////////////
-      // apply any boundary phase or twists
-      ////////////////////////////////////////////////////
       for (int mu = 0; mu < Nd; mu++) {
 
-	////////// boundary phase /////////////
+	      auto pha = Params.boundary_phases[mu];
-	auto pha = Params.boundary_phases[mu];
+	      scalar_type phase( real(pha),imag(pha) );
-	scalar_type phase( real(pha),imag(pha) );
 
-	int L   = GaugeGrid->GlobalDimensions()[mu];
+        int Lmu = GaugeGrid->GlobalDimensions()[mu] - 1;
-        int Lmu = L - 1;
 
         LatticeCoordinate(coor, mu);
 
         U = PeekIndex<LorentzIndex>(Umu, mu);
-
-	// apply any twists
-	RealD theta = Params.twist_n_2pi_L[mu] * 2*M_PI / L;
-	if ( theta != 0.0) { 
-	  scalar_type twphase(::cos(theta),::sin(theta));
-	  U = twphase*U;
-	  std::cout << GridLogMessage << " Twist ["<<mu<<"] "<< Params.twist_n_2pi_L[mu]<< " phase"<<phase <<std::endl;
-	}
-
         tmp = where(coor == Lmu, phase * U, U);
         PokeIndex<LorentzIndex>(Uds, tmp, mu);
 

Hadrons/A2AVectors.hpp

@@ -36,7 +36,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *                 Class to generate V & W all-to-all vectors                 *
+ *               Classes to generate V & W all-to-all vectors                 *
  ******************************************************************************/
 template <typename FImpl>
 class A2AVectorsSchurDiagTwo
@@ -70,42 +70,6 @@ private:
     SchurDiagTwoOperator<FMat, FermionField> op_;
 };
 
-/******************************************************************************
- *                  Methods for V & W all-to-all vectors I/O                  *
- ******************************************************************************/
-class A2AVectorsIo
-{
-public:
-    struct Record: Serializable
-    {
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Record,
-                                        unsigned int, index);
-        Record(void): index(0) {}
-    };
-public:
-    template <typename Field>
-    static void write(const std::string fileStem, std::vector<Field> &vec, 
-                      const bool multiFile, const int trajectory = -1);
-    template <typename Field>
-    static void read(std::vector<Field> &vec, const std::string fileStem,
-                     const bool multiFile, const int trajectory = -1);
-private:
-    static inline std::string vecFilename(const std::string stem, const int traj, 
-                                          const bool multiFile)
-    {
-        std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
-
-        if (multiFile)
-        {
-            return stem + t;
-        }
-        else
-        {
-            return stem + t + ".bin";
-        }
-    }
-};
-
 /******************************************************************************
  *               A2AVectorsSchurDiagTwo template implementation               *
  ******************************************************************************/
@@ -253,90 +217,6 @@ void A2AVectorsSchurDiagTwo<FImpl>::makeHighModeW5D(FermionField &wout_4d,
     }
 }
 
-/******************************************************************************
- *               all-to-all vectors I/O template implementation               *
- ******************************************************************************/
-template <typename Field>
-void A2AVectorsIo::write(const std::string fileStem, std::vector<Field> &vec, 
-                         const bool multiFile, const int trajectory)
-{
-    Record       record;
-    GridBase     *grid = vec[0]._grid;
-    ScidacWriter binWriter(grid->IsBoss());
-    std::string  filename = vecFilename(fileStem, trajectory, multiFile);
-
-    if (multiFile)
-    {
-        std::string fullFilename;
-
-        for (unsigned int i = 0; i < vec.size(); ++i)
-        {
-            fullFilename = filename + "/elem" + std::to_string(i) + ".bin";
-
-            LOG(Message) << "Writing vector " << i << std::endl;
-            makeFileDir(fullFilename, grid);
-            binWriter.open(fullFilename);
-            record.index = i;
-            binWriter.writeScidacFieldRecord(vec[i], record);
-            binWriter.close();
-        }
-    }
-    else
-    {
-        makeFileDir(filename, grid);
-        binWriter.open(filename);
-        for (unsigned int i = 0; i < vec.size(); ++i)
-        {
-            LOG(Message) << "Writing vector " << i << std::endl;
-            record.index = i;
-            binWriter.writeScidacFieldRecord(vec[i], record);
-        }
-        binWriter.close();
-    }
-}
-
-template <typename Field>
-void A2AVectorsIo::read(std::vector<Field> &vec, const std::string fileStem, 
-                        const bool multiFile, const int trajectory)
-{
-    Record       record;
-    ScidacReader binReader;
-    std::string  filename = vecFilename(fileStem, trajectory, multiFile);
-
-    if (multiFile)
-    {
-        std::string fullFilename;
-
-        for (unsigned int i = 0; i < vec.size(); ++i)
-        {
-            fullFilename = filename + "/elem" + std::to_string(i) + ".bin";
-
-            LOG(Message) << "Reading vector " << i << std::endl;
-            binReader.open(fullFilename);
-            binReader.readScidacFieldRecord(vec[i], record);
-            binReader.close();
-            if (record.index != i)
-            {
-                HADRONS_ERROR(Io, "vector index mismatch");
-            }
-        }
-    }
-    else
-    {
-        binReader.open(filename);
-        for (unsigned int i = 0; i < vec.size(); ++i)
-        {
-            LOG(Message) << "Reading vector " << i << std::endl;
-            binReader.readScidacFieldRecord(vec[i], record);
-            if (record.index != i)
-            {
-                HADRONS_ERROR(Io, "vector index mismatch");
-            }
-        }
-        binReader.close();
-    }
-}
-
 END_HADRONS_NAMESPACE
 
 #endif // A2A_Vectors_hpp_

Hadrons/DilutedNoise.hpp

@@ -7,7 +7,6 @@ Source file: Hadrons/DilutedNoise.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Vera Guelpers <Vera.Guelpers@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
@@ -77,22 +76,6 @@ private:
     unsigned int nt_;
 };
 
-template <typename FImpl>
-class FullVolumeSpinColorDiagonalNoise: public DilutedNoise<FImpl>
-{
-public:
-    typedef typename FImpl::FermionField FermionField;
-public:
-    // constructor/destructor
-    FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int n_src);
-    virtual ~FullVolumeSpinColorDiagonalNoise(void) = default;
-    // generate noise
-    virtual void generateNoise(GridParallelRNG &rng);
-private:
-    unsigned int nSrc_;
-};
-
-
 /******************************************************************************
  *                    DilutedNoise template implementation                    *
  ******************************************************************************/
@@ -203,47 +186,6 @@ void TimeDilutedSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rn
     }
 }
 
-/******************************************************************************
- *        FullVolumeSpinColorDiagonalNoise template implementation           *
- ******************************************************************************/
-template <typename FImpl>
-FullVolumeSpinColorDiagonalNoise<FImpl>::
-FullVolumeSpinColorDiagonalNoise(GridCartesian *g, unsigned int nSrc)
-: DilutedNoise<FImpl>(g, nSrc*Ns*FImpl::Dimension), nSrc_(nSrc)
-{}
-
-template <typename FImpl>
-void FullVolumeSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rng)
-{
-    typedef decltype(peekColour((*this)[0], 0)) SpinField;
-
-    auto                       &noise = *this;
-    auto                       g      = this->getGrid();
-    auto                       nd     = g->GlobalDimensions().size();
-    auto                       nc     = FImpl::Dimension;
-    Complex                    shift(1., 1.);
-    LatticeComplex             eta(g);
-    SpinField                  etas(g);
-    unsigned int               i = 0;
-
-    bernoulli(rng, eta);
-    eta = (2.*eta - shift)*(1./::sqrt(2.));
-    for (unsigned int n = 0; n < nSrc_; ++n)
-    {
-        for (unsigned int s = 0; s < Ns; ++s)
-        {
-            etas = zero;
-            pokeSpin(etas, eta, s);
-            for (unsigned int c = 0; c < nc; ++c)
-            {
-                noise[i] = zero;
-                pokeColour(noise[i], etas, c);
-                i++;
-            }
-        }
-    }
-}
-
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_DilutedNoise_hpp_

Hadrons/Modules.hpp

@@ -32,7 +32,6 @@
 #include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Hadrons/Modules/MGauge/StochEm.hpp>
 #include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
-#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
 #include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
 #include <Hadrons/Modules/MUtilities/RandomVectors.hpp>
 #include <Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
@@ -66,7 +65,6 @@
 #include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
 #include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
 #include <Hadrons/Modules/MIO/LoadNersc.hpp>
-#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
 #include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
 #include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
 #include <Hadrons/Modules/MIO/LoadBinary.hpp>

Hadrons/Modules/MAction/DWF.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TDWF<FIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TDWF<FIMPLF>;
-#endif

Hadrons/Modules/MAction/DWF.hpp

@@ -49,8 +49,7 @@ public:
                                     unsigned int, Ls,
                                     double      , mass,
                                     double      , M5,
-                                    std::string , boundary,
+                                    std::string , boundary);
-                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -74,9 +73,7 @@ protected:
 };
 
 MODULE_REGISTER_TMP(DWF, TDWF<FIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(DWFF, TDWF<FIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                        DWF template implementation                         *
@@ -120,9 +117,8 @@ void TDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    typename DomainWallFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
 }

Hadrons/Modules/MAction/MobiusDWF.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TMobiusDWF<FIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TMobiusDWF<FIMPLF>;
-#endif

Hadrons/Modules/MAction/MobiusDWF.hpp

@@ -49,8 +49,7 @@ public:
                                     double      , M5,
                                     double      , b,
                                     double      , c,
-                                    std::string , boundary,
+                                    std::string , boundary);
-                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -73,9 +72,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(MobiusDWF, TMobiusDWF<FIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(MobiusDWFF, TMobiusDWF<FIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                      TMobiusDWF implementation                             *
@@ -120,9 +117,8 @@ void TMobiusDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    typename MobiusFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                      implParams);

Hadrons/Modules/MAction/ScaledDWF.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TScaledDWF<FIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TScaledDWF<FIMPLF>;
-#endif

Hadrons/Modules/MAction/ScaledDWF.hpp

@@ -48,8 +48,7 @@ public:
                                     double      , mass,
                                     double      , M5,
                                     double      , scale,
-                                    std::string , boundary,
+                                    std::string , boundary);
-                                    std::string , twist);
 };
 
 template <typename FImpl>
@@ -72,9 +71,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(ScaledDWF, TScaledDWF<FIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(ScaledDWFF, TScaledDWF<FIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                      TScaledDWF implementation                             *
@@ -119,9 +116,8 @@ void TScaledDWF<FImpl>::setup(void)
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
-    typename ScaledShamirFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().scale,
                      implParams);

Hadrons/Modules/MAction/Wilson.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilson<FIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TWilson<FIMPLF>;
-#endif

Hadrons/Modules/MAction/Wilson.hpp

@@ -47,9 +47,7 @@ public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(WilsonPar,
                                     std::string, gauge,
                                     double     , mass,
-                                    std::string, boundary,
+                                    std::string, boundary);
-                                    std::string, string,
-                                    std::string, twist);
 };
 
 template <typename FImpl>
@@ -73,9 +71,7 @@ protected:
 };
 
 MODULE_REGISTER_TMP(Wilson, TWilson<FIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(WilsonF, TWilson<FIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                     TWilson template implementation                        *
@@ -115,9 +111,8 @@ void TWilson<FImpl>::setup(void)
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
-    typename WilsonFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
 }

Hadrons/Modules/MAction/WilsonClover.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilsonClover<FIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TWilsonClover<FIMPLF>;
-#endif

Hadrons/Modules/MAction/WilsonClover.hpp

@@ -51,8 +51,7 @@ public:
 				                    double     , csw_r,
 				                    double     , csw_t,
 				                    WilsonAnisotropyCoefficients ,clover_anisotropy,
-                                    std::string, boundary,
+                                    std::string, boundary
-                                    std::string, twist
 				    );
 };
 
@@ -76,9 +75,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(WilsonClover, TWilsonClover<FIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(WilsonCloverF, TWilsonClover<FIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                    TWilsonClover template implementation                   *
@@ -120,9 +117,8 @@ void TWilsonClover<FImpl>::setup(void)
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
-    typename WilsonCloverFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                      gridRb, par().mass, par().csw_r, par().csw_t, 
                      par().clover_anisotropy, implParams); 

Hadrons/Modules/MAction/ZMobiusDWF.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPL>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPLF>;
-#endif

Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -50,8 +50,7 @@ public:
                                     double                           , b,
                                     double                           , c,
                                     std::vector<std::complex<double>>, omega,
-                                    std::string                      , boundary,
+                                    std::string                      , boundary);
-                                    std::string                      , twist);
 };
 
 template <typename FImpl>
@@ -74,9 +73,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(ZMobiusDWFF, TZMobiusDWF<ZFIMPLF>, MAction);
-#endif
 
 /******************************************************************************
  *                     TZMobiusDWF implementation                             *
@@ -128,9 +125,8 @@ void TZMobiusDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
-    typename ZMobiusFermion<FImpl>::ImplParams implParams;
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    typename ZMobiusFermion<FImpl>::ImplParams implParams(boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
     envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, omega,
                      par().b, par().c, implParams);

Hadrons/Modules/MIO/LoadA2AVectors.cc

@@ -1,34 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MIO/LoadA2AVectors.cc
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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 <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MIO;
-
-template class Grid::Hadrons::MIO::TLoadA2AVectors<FIMPL>;

Hadrons/Modules/MIO/LoadA2AVectors.hpp

@@ -1,120 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MIO/LoadA2AVectors.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-
-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 Hadrons_MIO_LoadA2AVectors_hpp_
-#define Hadrons_MIO_LoadA2AVectors_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/A2AVectors.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                    Module to load all-to-all vectors                       *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MIO)
-
-class LoadA2AVectorsPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadA2AVectorsPar,
-                                    std::string,  filestem,
-                                    bool,         multiFile,
-                                    unsigned int, size);
-};
-
-template <typename FImpl>
-class TLoadA2AVectors: public Module<LoadA2AVectorsPar>
-{
-public:
-    FERM_TYPE_ALIASES(FImpl,);
-public:
-    // constructor
-    TLoadA2AVectors(const std::string name);
-    // destructor
-    virtual ~TLoadA2AVectors(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER_TMP(LoadA2AVectors, TLoadA2AVectors<FIMPL>, MIO);
-
-/******************************************************************************
- *                      TLoadA2AVectors implementation                        *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TLoadA2AVectors<FImpl>::TLoadA2AVectors(const std::string name)
-: Module<LoadA2AVectorsPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TLoadA2AVectors<FImpl>::getInput(void)
-{
-    std::vector<std::string> in;
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TLoadA2AVectors<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TLoadA2AVectors<FImpl>::setup(void)
-{
-    envCreate(std::vector<FermionField>, getName(), 1, par().size, 
-              envGetGrid(FermionField));
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TLoadA2AVectors<FImpl>::execute(void)
-{
-    auto &vec = envGet(std::vector<FermionField>, getName());
-
-    A2AVectorsIo::read(vec, par().filestem, par().multiFile, vm().getTrajectory());
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MIO_LoadA2AVectors_hpp_

Hadrons/Modules/MIO/LoadEigenPack.cc

@@ -32,6 +32,4 @@ using namespace Hadrons;
 using namespace MIO;
 
 template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>;
-#endif

Hadrons/Modules/MIO/LoadEigenPack.hpp

@@ -72,9 +72,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>), MIO);
-#endif
 
 /******************************************************************************
  *                    TLoadEigenPack implementation                           *

Hadrons/Modules/MNPR/Amputate.cc

@@ -6,7 +6,6 @@ Source file: Hadrons/Modules/MNPR/Amputate.cc
 
 Copyright (C) 2015-2018
 
-Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify

Hadrons/Modules/MNPR/Bilinear.cc

@@ -6,7 +6,6 @@ Source file: Hadrons/Modules/MNPR/Bilinear.cc
 
 Copyright (C) 2015-2018
 
-Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify

Hadrons/Modules/MNPR/FourQuark.cc

@@ -6,7 +6,6 @@ Source file: Hadrons/Modules/MNPR/FourQuark.cc
 
 Copyright (C) 2015-2018
 
-Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify

Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc

@@ -1,36 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Vera Guelpers <Vera.Guelpers@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 <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MNoise;
-
-template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal<FIMPL>;
-template class Grid::Hadrons::MNoise::TFullVolumeSpinColorDiagonal<ZFIMPL>;

Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp

@@ -1,121 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Vera Guelpers <Vera.Guelpers@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 */
-#ifndef Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_
-#define Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/DilutedNoise.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *             Generate full volume spin-color diagonal noise                *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MNoise)
-
-class FullVolumeSpinColorDiagonalPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(FullVolumeSpinColorDiagonalPar,
-                                    unsigned int, nsrc);
-};
-
-template <typename FImpl>
-class TFullVolumeSpinColorDiagonal: public Module<FullVolumeSpinColorDiagonalPar>
-{
-public:
-    FERM_TYPE_ALIASES(FImpl,);
-public:
-    // constructor
-    TFullVolumeSpinColorDiagonal(const std::string name);
-    // destructor
-    virtual ~TFullVolumeSpinColorDiagonal(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER_TMP(FullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal<FIMPL>, MNoise);
-MODULE_REGISTER_TMP(ZFullVolumeSpinColorDiagonal, TFullVolumeSpinColorDiagonal<ZFIMPL>, MNoise);
-
-/******************************************************************************
- *              TFullVolumeSpinColorDiagonal implementation                  *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TFullVolumeSpinColorDiagonal<FImpl>::TFullVolumeSpinColorDiagonal(const std::string name)
-: Module<FullVolumeSpinColorDiagonalPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TFullVolumeSpinColorDiagonal<FImpl>::getInput(void)
-{
-    std::vector<std::string> in;
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TFullVolumeSpinColorDiagonal<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TFullVolumeSpinColorDiagonal<FImpl>::setup(void)
-{
-    envCreateDerived(DilutedNoise<FImpl>, 
-                     FullVolumeSpinColorDiagonalNoise<FImpl>,
-                     getName(), 1, envGetGrid(FermionField), par().nsrc);
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TFullVolumeSpinColorDiagonal<FImpl>::execute(void)
-{
-    auto &noise = envGet(DilutedNoise<FImpl>, getName());
-    LOG(Message) << "Generating full volume, spin-color diagonal noise" << std::endl;
-    noise.generateNoise(rng4d());
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MNoise_FullVolumeSpinColorDiagonal_hpp_

Hadrons/Modules/MSolver/A2AVectors.hpp

@@ -51,9 +51,7 @@ public:
                                   std::string, noise,
                                   std::string, action,
                                   std::string, eigenPack,
-                                  std::string, solver,
+                                  std::string, solver);
-                                  std::string, output,
-                                  bool,        multiFile);
 };
 
 template <typename FImpl, typename Pack>
@@ -238,17 +236,6 @@ void TA2AVectors<FImpl, Pack>::execute(void)
         }
         stopTimer("W high mode");
     }
-
-    // I/O if necessary
-    if (!par().output.empty())
-    {
-        startTimer("V I/O");
-        A2AVectorsIo::write(par().output + "_v", v, par().multiFile, vm().getTrajectory());
-        stopTimer("V I/O");
-        startTimer("W I/O");
-        A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory());
-        stopTimer("W I/O");
-    }
 }
 
 END_MODULE_NAMESPACE

Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc

@@ -33,7 +33,4 @@ using namespace MSolver;
 
 template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
 template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<ZFIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
+
-template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS, FIMPLF>;
-template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<ZFIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS, ZFIMPLF>;
-#endif

Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -55,17 +55,17 @@ public:
                                     bool,          multiFile);
 };
 
-template <typename FImpl, int nBasis, typename FImplIo = FImpl>
+template <typename FImpl, int nBasis>
 class TLocalCoherenceLanczos: public Module<LocalCoherenceLanczosPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
     typedef LocalCoherenceLanczos<typename FImpl::SiteSpinor, 
                                   typename FImpl::SiteComplex, 
-                                  nBasis>                  LCL;
+                                  nBasis>                LCL;
-    typedef BaseFermionEigenPack<FImpl>                    BasePack;
+    typedef BaseFermionEigenPack<FImpl>                  BasePack;
-    typedef CoarseFermionEigenPack<FImpl, nBasis, FImplIo> CoarsePack;
+    typedef CoarseFermionEigenPack<FImpl, nBasis>        CoarsePack;
-    typedef HADRONS_DEFAULT_SCHUR_OP<FMat, FermionField>   SchurFMat;
+    typedef HADRONS_DEFAULT_SCHUR_OP<FMat, FermionField> SchurFMat;
 public:
     // constructor
     TLocalCoherenceLanczos(const std::string name);
@@ -82,31 +82,27 @@ public:
 
 MODULE_REGISTER_TMP(LocalCoherenceLanczos, ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
 MODULE_REGISTER_TMP(ZLocalCoherenceLanczos, ARG(TLocalCoherenceLanczos<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
-#ifdef GRID_DEFAULT_PRECISION_DOUBLE
-MODULE_REGISTER_TMP(LocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS, FIMPLF>), MSolver);
-MODULE_REGISTER_TMP(ZLocalCoherenceLanczosIo32, ARG(TLocalCoherenceLanczos<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS, ZFIMPLF>), MSolver);
-#endif
 
 /******************************************************************************
  *                 TLocalCoherenceLanczos implementation                      *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl, int nBasis, typename FImplIo>
+template <typename FImpl, int nBasis>
-TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::TLocalCoherenceLanczos(const std::string name)
+TLocalCoherenceLanczos<FImpl, nBasis>::TLocalCoherenceLanczos(const std::string name)
 : Module<LocalCoherenceLanczosPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl, int nBasis, typename FImplIo>
+template <typename FImpl, int nBasis>
-std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::getInput(void)
+std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getInput(void)
 {
     std::vector<std::string> in = {par().action};
     
     return in;
 }
 
-template <typename FImpl, int nBasis, typename FImplIo>
+template <typename FImpl, int nBasis>
-std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::getOutput(void)
+std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
@@ -114,8 +110,8 @@ std::vector<std::string> TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::getOutp
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl, int nBasis, typename FImplIo>
+template <typename FImpl, int nBasis>
-void TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::setup(void)
+void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
 {
     LOG(Message) << "Setting up local coherence Lanczos eigensolver for"
                  << " action '" << par().action << "' (" << nBasis
@@ -142,8 +138,8 @@ void TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::setup(void)
 }
 
 // execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl, int nBasis, typename FImplIo>
+template <typename FImpl, int nBasis>
-void TLocalCoherenceLanczos<FImpl, nBasis, FImplIo>::execute(void)
+void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
 {
     auto &finePar   = par().fineParams;
     auto &coarsePar = par().coarseParams;

Hadrons/Utilities/EigenPackCast.cc

@@ -35,7 +35,7 @@ using namespace Hadrons;
 template <typename FOut, typename FIn>
 void convert(const std::string outFilename, const std::string inFilename, 
              const unsigned int Ls, const bool rb, const unsigned int size, 
-             const bool multiFile, const bool testRead)
+             const bool multiFile)
 {
     assert(outFilename != inFilename);
     
@@ -102,7 +102,6 @@ void convert(const std::string outFilename, const std::string inFilename,
     LOG(Message) << "Out type      : " << typeName<FOut>() << std::endl;
     LOG(Message) << "#vectors      : " << size << std::endl;
     LOG(Message) << "Multifile     : " << (multiFile ? "yes" : "no") << std::endl;
-    LOG(Message) << "Test read     : " << (testRead ? "yes" : "no") << std::endl;
     if (multiFile)
     {
         for(unsigned int k = 0; k < size; ++k)
@@ -113,8 +112,6 @@ void convert(const std::string outFilename, const std::string inFilename,
             LOG(Message) << "==== Converting vector " << k << std::endl;
             LOG(Message) << "In : " << inV  << std::endl;
             LOG(Message) << "Out: " << outV << std::endl;
-            // conversion
-            LOG(Message) << "-- Doing conversion" << std::endl;
             makeFileDir(outV, gOut);
             binWriter.open(outV);
             binReader.open(inV);
@@ -124,20 +121,10 @@ void convert(const std::string outFilename, const std::string inFilename,
             EigenPackIo::writeElement<FIn, FOut>(binWriter, bufIn, eval, k, &bufOut, &testIn);
             binWriter.close();
             binReader.close();
-            // read test
-            if (testRead)
-            {
-                LOG(Message) << "-- Test read" << std::endl;
-                binReader.open(outV);
-                EigenPackIo::readElement<FOut>(bufOut, eval, k, binReader);
-                binReader.close();
-            }
         }
     }
     else
     {
-        // conversion
-        LOG(Message) << "-- Doing conversion" << std::endl;
         makeFileDir(outFilename, gOut);
         binWriter.open(outFilename);
         binReader.open(inFilename);
@@ -150,18 +137,6 @@ void convert(const std::string outFilename, const std::string inFilename,
         }
         binWriter.close();
         binReader.close();
-        // read test
-        if (testRead)
-        {
-            LOG(Message) << "-- Test read" << std::endl;
-            binReader.open(outFilename);
-            EigenPackIo::readHeader(record, binReader);
-            for(unsigned int k = 0; k < size; ++k)
-            {
-                EigenPackIo::readElement<FOut>(bufOut, eval, k, binReader);
-            }
-            binReader.close();
-        }
     }
 }
 
@@ -179,11 +154,11 @@ int main(int argc, char *argv[])
     // parse command line
     std::string  outFilename, inFilename;
     unsigned int size, Ls;
-    bool         rb, multiFile, testRead;
+    bool         rb, multiFile;
     
-    if (argc < 8)
+    if (argc < 7)
     {
-        std::cerr << "usage: " << argv[0] << " <out eigenpack> <in eigenpack> <Ls> <red-black {0|1}> <#vector> <multifile {0|1}> <test read {0|1}> [Grid options]";
+        std::cerr << "usage: " << argv[0] << " <out eigenpack> <in eigenpack> <Ls> <red-black (0|1)> <#vector> <multifile (0|1)> [Grid options]";
         std::cerr << std::endl;
         std::exit(EXIT_FAILURE);
     }
@@ -193,7 +168,6 @@ int main(int argc, char *argv[])
     rb          = (std::string(argv[4]) != "0");
     size        = std::stoi(std::string(argv[5]));
     multiFile   = (std::string(argv[6]) != "0");
-    testRead    = (std::string(argv[7]) != "0");
     
     // initialization
     Grid_init(&argc, &argv);
@@ -202,7 +176,7 @@ int main(int argc, char *argv[])
     // execution
     try
     {
-        convert<FOUT, FIN>(outFilename, inFilename, Ls, rb, size, multiFile, testRead);
+        convert<FOUT, FIN>(outFilename, inFilename, Ls, rb, size, multiFile);
     }
     catch (const std::exception& e)
     {

Hadrons/modules.inc

@@ -31,7 +31,6 @@ modules_cc =\
   Modules/MGauge/FundtoHirep.cc \
   Modules/MGauge/GaugeFix.cc \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
-  Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
   Modules/MUtilities/RandomVectors.cc \
   Modules/MUtilities/TestSeqGamma.cc \
   Modules/MUtilities/PrecisionCast.cc \
@@ -65,8 +64,7 @@ modules_cc =\
   Modules/MIO/LoadBinary.cc \
   Modules/MIO/LoadNersc.cc \
   Modules/MIO/LoadCoarseEigenPack.cc \
-  Modules/MIO/LoadCosmHol.cc \
+  Modules/MIO/LoadCosmHol.cc
-  Modules/MIO/LoadA2AVectors.cc
 
 modules_hpp =\
   Modules/MContraction/Baryon.hpp \
@@ -103,7 +101,6 @@ modules_hpp =\
   Modules/MGauge/FundtoHirep.hpp \
   Modules/MGauge/StochEm.hpp \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
-  Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
   Modules/MUtilities/PrecisionCast.hpp \
   Modules/MUtilities/RandomVectors.hpp \
   Modules/MUtilities/TestSeqGamma.hpp \
@@ -137,7 +134,6 @@ modules_hpp =\
   Modules/MScalarSUN/TrKinetic.hpp \
   Modules/MIO/LoadEigenPack.hpp \
   Modules/MIO/LoadNersc.hpp \
-  Modules/MIO/LoadA2AVectors.hpp \
   Modules/MIO/LoadCosmHol.hpp \
   Modules/MIO/LoadCoarseEigenPack.hpp \
   Modules/MIO/LoadBinary.hpp

benchmarks/Benchmark_IO_vs_dir.cc

@@ -3,6 +3,9 @@
 #define MSG std::cout << GridLogMessage
 #define SEP \
 "============================================================================="
+#ifndef BENCH_IO_LMAX
+#define BENCH_IO_LMAX 40
+#endif
 
 using namespace Grid;
 using namespace QCD;
@@ -38,7 +41,7 @@ int main (int argc, char ** argv)
   int64_t threads = GridThread::GetThreads();
   MSG << "Grid is setup to use " << threads << " threads" << std::endl;
   MSG << SEP << std::endl;
-  MSG << "Benchmark double precision Lime write" << std::endl;
+  MSG << "Benchmark Lime write" << std::endl;
   MSG << SEP << std::endl;
   for (auto &d: dir)
   {
@@ -46,8 +49,7 @@ int main (int argc, char ** argv)
     writeBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeWrite<LatticeFermion>, Ls, rb);
   }
 
-  MSG << SEP << std::endl;
+  MSG << "Benchmark Lime read" << std::endl;
-  MSG << "Benchmark double precision Lime read" << std::endl;
   MSG << SEP << std::endl;
   for (auto &d: dir)
   {
@@ -55,24 +57,6 @@ int main (int argc, char ** argv)
     readBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeRead<LatticeFermion>, Ls, rb);
   }
 
-  MSG << SEP << std::endl;
-  MSG << "Benchmark single precision Lime write" << std::endl;
-  MSG << SEP << std::endl;
-  for (auto &d: dir)
-  {
-    MSG << "-- Directory " << d << std::endl;
-    writeBenchmark<LatticeFermionF>(GridDefaultLatt(), d + "/ioBench", limeWrite<LatticeFermionF>, Ls, rb);
-  }
-
-  MSG << SEP << std::endl;
-  MSG << "Benchmark single precision Lime read" << std::endl;
-  MSG << SEP << std::endl;
-  for (auto &d: dir)
-  {
-    MSG << "-- Directory " << d << std::endl;
-    readBenchmark<LatticeFermionF>(GridDefaultLatt(), d + "/ioBench", limeRead<LatticeFermionF>, Ls, rb);
-  }
-
   Grid_finalize();
 
   return EXIT_SUCCESS;

configure.ac

@@ -485,7 +485,6 @@ DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg])
 
 ############### Ouput
 cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
-GRID_CXX="$CXX"
 GRID_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
 GRID_LDFLAGS="$AM_LDFLAGS $LDFLAGS"
 GRID_LIBS=$LIBS
@@ -498,7 +497,6 @@ AM_LDFLAGS="-L${cwd}/Grid $AM_LDFLAGS"
 AC_SUBST([AM_CFLAGS])
 AC_SUBST([AM_CXXFLAGS])
 AC_SUBST([AM_LDFLAGS])
-AC_SUBST([GRID_CXX])
 AC_SUBST([GRID_CXXFLAGS])
 AC_SUBST([GRID_LDFLAGS])
 AC_SUBST([GRID_LIBS])

grid-config.in

@@ -61,10 +61,6 @@ while test $# -gt 0; do
       echo @GRID_CXXFLAGS@
     ;;
     
-    --cxx)
-      echo @GRID_CXX@
-    ;;
-    
     --ldflags)
       echo @GRID_LDFLAGS@
     ;;

tests/debug/Test_split_laplacian.cc

@@ -1,104 +0,0 @@
-    /*************************************************************************************
-
-    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>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-
-int main (int argc, char ** argv)
-{
-  typedef LatticeComplex ComplexField; 
-
-  Grid_init(&argc,&argv);
-
-  std::vector<int> latt_size   = GridDefaultLatt();
-  int nd   = latt_size.size();
-  int ndm1 = nd-1;
-
-  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::cout << " Full " << GridCmdVectorIntToString(latt_size)  << " subgrid"         <<std::endl;
-  std::cout << " Full " << GridCmdVectorIntToString(mpi_layout) << " sub communicator"<<std::endl;
-  std::cout << " Full " << GridCmdVectorIntToString(simd_layout)<< " simd layout "    <<std::endl;
-
-  GridCartesian         * GridN = new GridCartesian(latt_size,
-						    simd_layout,
-						    mpi_layout);
-
-  std::vector<int> latt_m  = latt_size;   latt_m[nd-1] = 1;
-  std::vector<int> mpi_m   = mpi_layout;  mpi_m [nd-1] = 1;
-  std::vector<int> simd_m  = GridDefaultSimd(ndm1,vComplex::Nsimd()); simd_m.push_back(1);
-
-
-  std::cout << " Requesting " << GridCmdVectorIntToString(latt_m)<< " subgrid"         <<std::endl;
-  std::cout << " Requesting " << GridCmdVectorIntToString(mpi_m) << " sub communicator"<<std::endl;
-  std::cout << " Requesting " << GridCmdVectorIntToString(simd_m)<< " simd layout "    <<std::endl;
-  GridCartesian         * Grid_m = new GridCartesian(latt_m,
-						     simd_m,
-						     mpi_m,
-						     *GridN); 
-
-  Complex C(1.0);
-  Complex tmp;
-
-  ComplexField Full(GridN); Full = C;
-  ComplexField Full_cpy(GridN);
-  ComplexField Split(Grid_m);Split= C;
-
-  std::cout << GridLogMessage<< " Full  volume "<< norm2(Full) <<std::endl;
-  std::cout << GridLogMessage<< " Split volume "<< norm2(Split) <<std::endl;
-
-  tmp=C;
-  GridN->GlobalSum(tmp);
-  std::cout << GridLogMessage<< " Full  nodes "<< tmp <<std::endl;
-
-  tmp=C;
-  Grid_m->GlobalSum(tmp);
-  std::cout << GridLogMessage<< " Split nodes "<< tmp <<std::endl;
-  GridN->Barrier();
-
-  auto local_latt = GridN->LocalDimensions();
-
-  Full_cpy = zero;
-  std::vector<int> seeds({1,2,3,4});
-  GridParallelRNG          RNG(GridN);  RNG.SeedFixedIntegers(seeds);
-
-  random(RNG,Full);
-  for(int t=0;t<local_latt[nd-1];t++){
-    ExtractSliceLocal(Split,Full,0,t,Tp);
-    InsertSliceLocal (Split,Full_cpy,0,t,Tp);
-  }
-  Full_cpy = Full_cpy - Full;
-  std::cout << " NormFull " << norm2(Full)<<std::endl;
-  std::cout << " NormDiff " << norm2(Full_cpy)<<std::endl;
-  Grid_finalize();
-}

tests/hadrons/Test_QED.cc

@@ -72,7 +72,6 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
-    std::string twist    = "0. 0. 0. 0.";
 
     //stochastic photon field
     MGauge::StochEm::Par photonPar;
@@ -91,7 +90,6 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
-        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
 
         

tests/hadrons/Test_hadrons.hpp

@@ -126,7 +126,6 @@ inline void makeWilsonAction(Application &application, std::string actionName,
         actionPar.gauge = gaugeField;
         actionPar.mass  = mass;
         actionPar.boundary = boundary;
-        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::Wilson>(actionName, actionPar);
     }
 }
@@ -155,7 +154,6 @@ inline void makeDWFAction(Application &application, std::string actionName,
         actionPar.M5    = M5;
         actionPar.mass  = mass;
         actionPar.boundary = boundary;
-        actionPar.twist = "0. 0. 0. 0.";
         application.createModule<MAction::DWF>(actionName, actionPar);
     }
 }

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -66,7 +66,6 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
-    std::string twist = "0. 0. 0. 0.";
 
     // sink
     MSink::Point::Par sinkPar;
@@ -81,7 +80,6 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
-        actionPar.twist = twist;
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
         
         // solvers

tests/hadrons/Test_hadrons_spectrum.cc

@@ -72,7 +72,6 @@ int main(int argc, char *argv[])
     
     // set fermion boundary conditions to be periodic space, antiperiodic time.
     std::string boundary = "1 1 1 -1";
-    std::string twist = "0. 0. 0. 0.";
 
     for (unsigned int i = 0; i < flavour.size(); ++i)
     {
@@ -83,7 +82,6 @@ int main(int argc, char *argv[])
         actionPar.M5    = 1.8;
         actionPar.mass  = mass[i];
         actionPar.boundary = boundary;
-        actionPar.twist = twist;
         application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
         
         // solvers

tests/solver/Test_dwf_mrhs_cg_mpi.cc

@@ -38,7 +38,6 @@ int main (int argc, char ** argv)
   typedef typename DomainWallFermionR::ComplexField ComplexField; 
   typename DomainWallFermionR::ImplParams params; 
 
-  double stp=1.0e-5;
   const int Ls=4;
 
   Grid_init(&argc,&argv);
@@ -198,7 +197,7 @@ int main (int argc, char ** argv)
 
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
   MdagMLinearOperator<DomainWallFermionR,FermionField> HermOpCk(Dchk);
-  ConjugateGradient<FermionField> CG((stp),10000);
+  ConjugateGradient<FermionField> CG((1.0e-2),10000);
   s_res = zero;
   CG(HermOp,s_src,s_res);
 
@@ -228,11 +227,5 @@ 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

@@ -1,220 +0,0 @@
-   /*************************************************************************************
-
-    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-8;
-  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;
-  GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
-  for(int s=0;s<nrhs;s++) {
-    random(pRNG5,src[s]);
-    std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
-  }
-
-  std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
-
-  LatticeGaugeField Umu(UGrid); 
-
-  if(0) { 
-    FieldMetaData header;
-    std::string file("./lat.in");
-    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;
-  } 
-
-  /////////////////
-  // 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;
-
-  ///////////////////////////////////////////////////////////////
-  // Set up N-solvers as trivially parallel
-  ///////////////////////////////////////////////////////////////
-  std::cout << GridLogMessage << " Building the solvers"<<std::endl;
-  //  RealD mass=0.00107;
-  RealD mass=0.1;
-  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;
-
-  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<<"]  "<< std::sqrt(norm2(tmp)/norm2(src[n]))<<std::endl;
-  }
-
-
-  for(int s=0;s<nrhs;s++){
-    result[s]=zero;
-  }
-
-  /////////////////////////////////////////////////////////////
-  // Try block CG
-  /////////////////////////////////////////////////////////////
-  int blockDim = 0;//not used for BlockCGVec
-  BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
-  {
-    BCGV(HermOpCk,src,result);
-  }
-
-  
-  
-  Grid_finalize();
-}

tests/solver/Test_mobius_bcg_nosplit.cc

@@ -1,144 +0,0 @@
-   /*************************************************************************************
-
-    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 DomainWallFermionR::FermionField FermionField; 
-  typedef typename DomainWallFermionR::ComplexField ComplexField; 
-  typename DomainWallFermionR::ImplParams params; 
-
-  const int Ls=16;
-
-  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<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);
- 
-  double stp = 1.e-8;
-  int nrhs = 2;
-
-  ///////////////////////////////////////////////
-  // 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;
-  GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
-  for(int s=0;s<nrhs;s++) {
-    random(pRNG5,src[s]);
-    std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
-  }
-
-  std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
-
-  LatticeGaugeField Umu(UGrid); 
-
-  int conf = 0;
-  if(conf==0) { 
-    FieldMetaData header;
-    std::string file("./lat.in");
-    NerscIO::readConfiguration(Umu,header,file);
-    std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
-  } else if (conf==1){
-    GridParallelRNG pRNG(UGrid );  
-
-    pRNG.SeedFixedIntegers(seeds);
-    SU3::HotConfiguration(pRNG,Umu);
-    std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
-  } else {
-    SU3::ColdConfiguration(Umu);
-    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
-  }
-
-  ///////////////////////////////////////////////////////////////
-  // Set up N-solvers as trivially parallel
-  ///////////////////////////////////////////////////////////////
-  std::cout << GridLogMessage << " Building the solvers"<<std::endl;
-  RealD mass=0.01;
-  RealD M5=1.8;
-  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
-
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-
-  MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
-  ConjugateGradient<FermionField> CG((stp),100000);
-
-  for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
-    CG(HermOp,src[rhs],result[rhs]);
-  }
-
-  for(int rhs=0;rhs<1;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  /////////////////////////////////////////////////////////////
-  // Try block CG
-  /////////////////////////////////////////////////////////////
-  int blockDim = 0;//not used for BlockCGVec
-  for(int s=0;s<nrhs;s++){
-      result[s]=zero;
-   }
-  BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
-  {
-    BCGV(HermOp,src,result);
-  }
-  
-  for(int rhs=0;rhs<nrhs;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  Grid_finalize();
-}

tests/solver/Test_mobius_bcg_phys_nosplit.cc

@@ -1,148 +0,0 @@
-   /*************************************************************************************
-
-    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 DomainWallFermionR::FermionField FermionField; 
-  typedef typename DomainWallFermionR::ComplexField ComplexField; 
-  typename DomainWallFermionR::ImplParams params; 
-
-  const int Ls=16;
-
-  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<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);
- 
-  double stp = 1.e-8;
-  int nrhs = 2;
-
-  ///////////////////////////////////////////////
-  // Set up the problem as a 4d spreadout job
-  ///////////////////////////////////////////////
-  std::vector<int> seeds({1,2,3,4});
-
-  std::vector<FermionField> src4(nrhs,UGrid);
-  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;
-  GridParallelRNG pRNG4(UGrid);  pRNG4.SeedFixedIntegers(seeds);
-  for(int s=0;s<nrhs;s++) {
-    random(pRNG4,src4[s]);
-    std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
-  }
-
-  std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
-
-  LatticeGaugeField Umu(UGrid); 
-
-  int conf = 0;
-  if(conf==0) { 
-    FieldMetaData header;
-    std::string file("./lat.in");
-    NerscIO::readConfiguration(Umu,header,file);
-    std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
-  } else if (conf==1){
-    GridParallelRNG pRNG(UGrid );  
-
-    pRNG.SeedFixedIntegers(seeds);
-    SU3::HotConfiguration(pRNG,Umu);
-    std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
-  } else {
-    SU3::ColdConfiguration(Umu);
-    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
-  }
-
-  ///////////////////////////////////////////////////////////////
-  // Set up N-solvers as trivially parallel
-  ///////////////////////////////////////////////////////////////
-  std::cout << GridLogMessage << " Building the solvers"<<std::endl;
-  RealD mass=0.01;
-  RealD M5=1.8;
-  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
-  for(int s=0;s<nrhs;s++) {
-    Ddwf.ImportPhysicalFermionSource(src4[s],src[s]);
-  }
-
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-
-  MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
-  ConjugateGradient<FermionField> CG((stp),100000);
-
-  for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
-    //    CG(HermOp,src[rhs],result[rhs]);
-  }
-
-  for(int rhs=0;rhs<1;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  /////////////////////////////////////////////////////////////
-  // Try block CG
-  /////////////////////////////////////////////////////////////
-  int blockDim = 0;//not used for BlockCGVec
-  for(int s=0;s<nrhs;s++){
-    result[s]=zero;
-  }
-  BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
-  {
-    BCGV(HermOp,src,result);
-  }
-  
-  for(int rhs=0;rhs<nrhs;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  Grid_finalize();
-}

tests/solver/Test_mobius_bcg_prec_nosplit.cc

@@ -1,147 +0,0 @@
-   /*************************************************************************************
-
-    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 DomainWallFermionR::FermionField FermionField; 
-  typedef typename DomainWallFermionR::ComplexField ComplexField; 
-  typename DomainWallFermionR::ImplParams params; 
-
-  const int Ls=16;
-
-  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<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);
- 
-  double stp = 1.e-8;
-  int nrhs = 4;
-
-  ///////////////////////////////////////////////
-  // 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;
-  GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
-  for(int s=0;s<nrhs;s++) {
-    random(pRNG5,src[s]);
-    std::cout << GridLogMessage << " src ["<<s<<"] "<<norm2(src[s])<<std::endl;
-  }
-
-  std::cout << GridLogMessage << "Intialised the Fermion Fields"<<std::endl;
-
-  LatticeGaugeField Umu(UGrid); 
-
-  int conf = 2;
-  if(conf==0) { 
-    FieldMetaData header;
-    std::string file("./lat.in");
-    NerscIO::readConfiguration(Umu,header,file);
-    std::cout << GridLogMessage << " Config "<<file<<" successfully read" <<std::endl;
-  } else if (conf==1){
-    GridParallelRNG pRNG(UGrid );  
-
-    pRNG.SeedFixedIntegers(seeds);
-    SU3::HotConfiguration(pRNG,Umu);
-    std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
-  } else {
-    SU3::ColdConfiguration(Umu);
-    std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
-  }
-
-  ///////////////////////////////////////////////////////////////
-  // Set up N-solvers as trivially parallel
-  ///////////////////////////////////////////////////////////////
-  std::cout << GridLogMessage << " Building the solvers"<<std::endl;
-  RealD mass=0.01;
-  RealD M5=1.8;
-  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*rbGrid,mass,M5,params);
-
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::cout << GridLogMessage << " Calling DWF CG "<<std::endl;
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-
-  MdagMLinearOperator<DomainWallFermionR,FermionField> HermOp(Ddwf);
-  ConjugateGradient<FermionField> CG((stp),100000);
-
-  for(int rhs=0;rhs<1;rhs++){
-    result[rhs] = zero;
-    CG(HermOp,src[rhs],result[rhs]);
-  }
-
-  for(int rhs=0;rhs<1;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  /////////////////////////////////////////////////////////////
-  // Try block CG
-  /////////////////////////////////////////////////////////////
-  int blockDim = 0;//not used for BlockCGVec
-  for(int s=0;s<nrhs;s++){
-    result[s]=zero;
-  }
-
-
-  {
-    BlockConjugateGradient<FermionField>    BCGV  (BlockCGrQVec,blockDim,stp,100000);
-    SchurRedBlackDiagTwoSolve<FermionField> SchurSolver(BCGV);
-    SchurSolver(Ddwf,src,result);
-  }
-  
-  for(int rhs=0;rhs<nrhs;rhs++){
-    std::cout << " Result["<<rhs<<"] norm = "<<norm2(result[rhs])<<std::endl;
-  }
-
-  Grid_finalize();
-}

tests/solver/Test_staggered_block_cg_prec.cc

@@ -67,22 +67,7 @@ int main (int argc, char ** argv)
   GridParallelRNG pRNG(UGrid );  pRNG.SeedFixedIntegers(seeds);
   GridParallelRNG pRNG5(FGrid);  pRNG5.SeedFixedIntegers(seeds);
 
-  FermionField src(FGrid);
+  FermionField src(FGrid); random(pRNG5,src);
-  FermionField tt(FGrid);
-#if 1
-  random(pRNG5,src);
-#else
-  src=zero;
-  ComplexField coor(FGrid);
-  LatticeCoordinate(coor,0);
-  for(int ss=0;ss<FGrid->oSites();ss++){
-    src._odata[ss]()()(0)=coor._odata[ss]()()();
-  }
-  LatticeCoordinate(coor,1);
-  for(int ss=0;ss<FGrid->oSites();ss++){
-    src._odata[ss]()()(0)+=coor._odata[ss]()()();
-  }
-#endif
   FermionField src_o(FrbGrid);   pickCheckerboard(Odd,src_o,src);
   FermionField result_o(FrbGrid); result_o=zero; 
   RealD nrm = norm2(src);
@@ -104,8 +89,7 @@ int main (int argc, char ** argv)
   ConjugateGradient<FermionField> CG(1.0e-8,10000);
   int blockDim = 0;
   BlockConjugateGradient<FermionField>    BCGrQ(BlockCGrQ,blockDim,1.0e-8,10000);
-  BlockConjugateGradient<FermionField>    BCG  (BlockCGrQ,blockDim,1.0e-8,10000);
+  BlockConjugateGradient<FermionField>    BCG  (BlockCG,blockDim,1.0e-8,10000);
-  BlockConjugateGradient<FermionField>    BCGv (BlockCGrQVec,blockDim,1.0e-8,10000);
   BlockConjugateGradient<FermionField>    mCG  (CGmultiRHS,blockDim,1.0e-8,10000);
 
   std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
@@ -174,7 +158,7 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
 
   std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
-  std::cout << GridLogMessage << " Calling Block CGrQ for "<<Ls <<" right hand sides" <<std::endl;
+  std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
   std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
   Ds.ZeroCounters();
   result_o=zero;
@@ -192,49 +176,6 @@ int main (int argc, char ** argv)
   Ds.Report();
   std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
 
-  std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
-  std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
-  std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
-  Ds.ZeroCounters();
-  result_o=zero;
-  {
-    double t1=usecond();
-    BCG(HermOp,src_o,result_o);
-    double t2=usecond();
-    double ncall=BCGrQ.IterationsToComplete*Ls;
-    double flops = deodoe_flops * ncall;
-    std::cout<<GridLogMessage << "usec    =   "<< (t2-t1)<<std::endl;
-    std::cout<<GridLogMessage << "flops   =   "<< flops<<std::endl;
-    std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t2-t1)<<std::endl;
-    HermOp.Report();
-  }
-  Ds.Report();
-  std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
-
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::cout << GridLogMessage << " Calling BCGvec "<<std::endl;
-  std::cout << GridLogMessage << "****************************************************************** "<<std::endl;
-  std::vector<FermionField> src_v   (Ls,UrbGrid);
-  std::vector<FermionField> result_v(Ls,UrbGrid);
-  for(int s=0;s<Ls;s++) result_v[s] = zero;
-  for(int s=0;s<Ls;s++) {
-    FermionField src4(UGrid);
-    ExtractSlice(src4,src,s,0);
-    pickCheckerboard(Odd,src_v[s],src4);  
-  }
-
-  {
-    double t1=usecond();
-    BCGv(HermOp4d,src_v,result_v);
-    double t2=usecond();
-    double ncall=BCGv.IterationsToComplete*Ls;
-    double flops = deodoe_flops * ncall;
-    std::cout<<GridLogMessage << "usec    =   "<< (t2-t1)<<std::endl;
-    std::cout<<GridLogMessage << "flops   =   "<< flops<<std::endl;
-    std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t2-t1)<<std::endl;
-    //    HermOp4d.Report();
-  }
-
 
   Grid_finalize();
 }