Merge 32e6d58356 into aa67a5b095

BLAS_benchmark/BatchBlasBench.cc

@@ -12,15 +12,13 @@
 #include <iostream>
 #include <sys/time.h>
 
-#define GRID_SYCL
-#undef  GRID_HIP
-#undef  GRID_CUDA
 
 #ifdef GRID_HIP
 #include <hipblas/hipblas.h>
 #endif
 #ifdef GRID_CUDA
 #include <cublas_v2.h>
+
 #endif
 #ifdef GRID_SYCL
 #include <oneapi/mkl.hpp>
@@ -45,6 +43,90 @@ inline void acceleratorFreeDevice(void *ptr,size_t bytes){free(ptr,*theAccelerat
 inline void acceleratorMemSet(void *base,int value,size_t bytes) { theAccelerator->memset(base,value,bytes); theAccelerator->wait();}
 inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { theAccelerator->memcpy(to,from,bytes); theAccelerator->wait();}
 inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ theAccelerator->memcpy(to,from,bytes); theAccelerator->wait();}
+#define accelerator_barrier(dummy) { theAccelerator->wait(); }
+#endif
+
+#ifdef GRID_HIP
+hipStream_t copyStream;
+hipStream_t computeStream;
+void acceleratorInit(void)
+{
+  int device = 0;
+  auto discard = hipSetDevice(device);
+  discard = hipStreamCreate(&copyStream);
+  discard = hipStreamCreate(&computeStream);
+  printf("AcceleratorHIPInit\n");
+}
+inline void *acceleratorAllocDevice(size_t bytes)
+{
+  void *ptr=NULL;
+  auto err = hipMalloc((void **)&ptr,bytes);
+  if( err != hipSuccess ) {
+    ptr = (void *) NULL;
+    fprintf(stderr," hipMalloc failed for %ld %s \n",bytes,hipGetErrorString(err)); fflush(stderr);
+  }
+  return ptr;
+};
+inline void acceleratorFreeDevice(void *ptr,size_t bytes){ auto discard=hipFree(ptr);};
+inline void acceleratorFreeDevice(void *ptr){ auto discard=hipFree(ptr);};
+inline void acceleratorMemSet(void *base,int value,size_t bytes) { auto discard=hipMemset(base,value,bytes);}
+inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { auto discard=hipMemcpy(to,from,bytes, hipMemcpyHostToDevice);}
+inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ auto discard=hipMemcpy(to,from,bytes, hipMemcpyDeviceToHost);}
+#define accelerator_barrier(dummy)				\
+  {								\
+    auto tmp=hipStreamSynchronize(computeStream);		\
+    auto err = hipGetLastError();				\
+    if ( err != hipSuccess ) {					\
+      printf("After hipDeviceSynchronize() : HIP error %s \n", hipGetErrorString( err )); \
+      puts(__FILE__);							\
+      printf("Line %d\n",__LINE__);				\
+      exit(0);							\
+    }								\
+  }
+
+#endif
+
+#ifdef GRID_CUDA
+cudaStream_t copyStream;
+cudaStream_t computeStream;
+void acceleratorInit(void)
+{
+  int device = 0;
+  cudaSetDevice(device);
+  cudaStreamCreate(&copyStream);
+  cudaStreamCreate(&computeStream);
+}
+inline void *acceleratorAllocDevice(size_t bytes)
+{
+  void *ptr=NULL;
+  auto err = cudaMalloc((void **)&ptr,bytes);
+  if( err != cudaSuccess ) {
+    ptr = (void *) NULL;
+    printf(" cudaMalloc failed for %d %s \n",bytes,cudaGetErrorString(err));
+  }
+  return ptr;
+};
+inline void acceleratorFreeShared(void *ptr){ cudaFree(ptr);};
+inline void acceleratorFreeDevice(void *ptr){ cudaFree(ptr);};
+inline void acceleratorCopyToDevice(void *from,void *to,size_t bytes)  { cudaMemcpy(to,from,bytes, cudaMemcpyHostToDevice);}
+inline void acceleratorCopyFromDevice(void *from,void *to,size_t bytes){ cudaMemcpy(to,from,bytes, cudaMemcpyDeviceToHost);}
+inline void acceleratorMemSet(void *base,int value,size_t bytes) { cudaMemset(base,value,bytes);}
+#define accelerator_barrier(dummy)					\
+  {									\
+    cudaStreamSynchronize(computeStream);				\
+    cudaError err = cudaGetLastError();					\
+    if ( cudaSuccess != err ) {						\
+      printf("accelerator_barrier(): Cuda error %s \n",			\
+	     cudaGetErrorString( err ));				\
+      printf("File %s Line %d\n",__FILE__,__LINE__);			\
+      fflush(stdout);							\
+      if (acceleratorAbortOnGpuError) assert(err==cudaSuccess);		\
+    }									\
+  }
+
+#endif
+
+
 template<class T> void acceleratorPut(T& dev,T&host)
 {
   acceleratorCopyToDevice(&host,&dev,sizeof(T));
@@ -55,9 +137,6 @@ template<class T> T acceleratorGet(T& dev)
   acceleratorCopyFromDevice(&dev,&host,sizeof(T));
   return host;
 }
-#define accelerator_barrier(dummy) { theAccelerator->wait(); }
-#endif
-
 
 /**************************************************************
  * Allocator
@@ -210,7 +289,270 @@ public:
     gridblasHandle->wait();
 #endif
   }
+
+
+  /////////////////////////////////////////////////////////////
+  // Single matrix GEMM -- fp64 and fp32
+  /////////////////////////////////////////////////////////////
+  void gemm(GridBLASOperation_t OpA,
+	    GridBLASOperation_t OpB,
+	    int m,int n, int k,
+	    ComplexD alpha,
+	    ComplexD* Amk,  // Device pointer
+	    ComplexD* Bkn,
+	    ComplexD beta,
+	    ComplexD* Cmn)
+  {
+    RealD t2=usecond();
+
+    assert(OpA!=GridBLAS_OP_T); // Complex case expect no transpose
+    assert(OpB!=GridBLAS_OP_T);
+
+    int lda = m; // m x k column major
+    int ldb = k; // k x n column major
+    int ldc = m; // m x b column major
+    if(OpA!=GridBLAS_OP_N)
+      lda = k;
+    if(OpB!=GridBLAS_OP_N)
+      ldb = n;
+    
+    static deviceVector<ComplexD> alpha_p(1);
+    static deviceVector<ComplexD> beta_p(1);
+    // can prestore the 1 and the zero on device
+    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(ComplexD));
+    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(ComplexD));
+    RealD t0=usecond();
+
+#ifdef GRID_HIP
+    hipblasOperation_t hOpA;
+    hipblasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
+    auto err = hipblasZgemm(gridblasHandle,
+			    hOpA,
+			    hOpB,
+			    m,n,k,
+			    (hipblasDoubleComplex *) &alpha_p[0],
+			    (hipblasDoubleComplex *) Amk, lda,
+			    (hipblasDoubleComplex *) Bkn, ldb,
+			    (hipblasDoubleComplex *) &beta_p[0],
+			    (hipblasDoubleComplex *) Cmn, ldc);
+    assert(err==HIPBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_CUDA
+    cublasOperation_t hOpA;
+    cublasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
+    auto err = cublasZgemm(gridblasHandle,
+			   hOpA,
+			   hOpB,
+			   m,n,k,
+			   (cuDoubleComplex *) &alpha_p[0],
+			   (cuDoubleComplex *) Amk, lda,
+			   (cuDoubleComplex *) Bkn, ldb,
+			   (cuDoubleComplex *) &beta_p[0],
+			   (cuDoubleComplex *) Cmn, ldc);
+    assert(err==CUBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_SYCL
+      int64_t m64=m;
+      int64_t n64=n;
+      int64_t k64=k;
+      int64_t lda64=lda;
+      int64_t ldb64=ldb;
+      int64_t ldc64=ldc;
+
+      oneapi::mkl::transpose iOpA;
+      oneapi::mkl::transpose iOpB;
+      
+      if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
+      if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
+      if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
+      if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
+      if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
+      if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
+
+      oneapi::mkl::blas::column_major::gemm(*gridblasHandle,
+					    iOpA,
+					    iOpB,
+					    m64,n64,k64,
+					    (ComplexD *) &alpha_p[0],
+					    (const ComplexD *)Amk, (int64_t )lda64,
+					    (const ComplexD *)Bkn, (int64_t )ldb64,
+					    (ComplexD *) &beta_p[0],
+					    (ComplexD *)Cmn, (int64_t)ldc64);
+      synchronise();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    // Need a default/reference implementation; use Eigen
+      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+	Eigen::Map<Eigen::MatrixXcd> eAmk(Amk,m,k);
+	Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn,k,n);
+	Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+	Eigen::Map<Eigen::MatrixXcd> eAmk(Amk,k,m);
+	Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn,k,n);
+	Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn ;
+      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
+	Eigen::Map<Eigen::MatrixXcd> eAmk(Amk,m,k);
+	Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn,n,k);
+	Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint() ;
+      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
+	Eigen::Map<Eigen::MatrixXcd> eAmk(Amk,k,m);
+	Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn,n,k);
+	Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint() ;
+      } else { 
+	assert(0);
+      }
+#endif
+     RealD t1=usecond();
+     RealD flops = 8.0*m*n*k;
+     RealD bytes = 1.0*sizeof(ComplexD)*(m*k+k*n+m*n);
+  }
+  void gemm(GridBLASOperation_t OpA,
+	    GridBLASOperation_t OpB,
+	    int m,int n, int k,
+	    ComplexF alpha,
+	    ComplexF* Amk,  // Device pointer
+	    ComplexF* Bkn,
+	    ComplexF beta,
+	    ComplexF* Cmn)
+  {
+    RealD t2=usecond();
+
+    assert(OpA!=GridBLAS_OP_T); // Complex case expect no transpose
+    assert(OpB!=GridBLAS_OP_T);
+
+    int lda = m; // m x k column major
+    int ldb = k; // k x n column major
+    int ldc = m; // m x b column major
+    if(OpA!=GridBLAS_OP_N)
+      lda = k;
+    if(OpB!=GridBLAS_OP_N)
+      ldb = n;
+    
+    static deviceVector<ComplexF> alpha_p(1);
+    static deviceVector<ComplexF> beta_p(1);
+    // can prestore the 1 and the zero on device
+    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(ComplexF));
+    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(ComplexF));
+    RealD t0=usecond();
+
+#ifdef GRID_HIP
+    hipblasOperation_t hOpA;
+    hipblasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
+    auto err = hipblasCgemm(gridblasHandle,
+			    hOpA,
+			    hOpB,
+			    m,n,k,
+			    (hipblasComplex *) &alpha_p[0],
+			    (hipblasComplex *) Amk, lda,
+			    (hipblasComplex *) Bkn, ldb,
+			    (hipblasComplex *) &beta_p[0],
+			    (hipblasComplex *) Cmn, ldc);
+    assert(err==HIPBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_CUDA
+    cublasOperation_t hOpA;
+    cublasOperation_t hOpB;
+    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
+    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
+    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
+    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
+    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
+    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
+    auto err = cublasCgemm(gridblasHandle,
+			   hOpA,
+			   hOpB,
+			   m,n,k,
+			   (cuComplex *) &alpha_p[0],
+			   (cuComplex *) Amk, lda,
+			   (cuComplex *) Bkn, ldb,
+			   (cuComplex *) &beta_p[0],
+			   (cuComplex *) Cmn, ldc);
+    assert(err==CUBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_SYCL
+      int64_t m64=m;
+      int64_t n64=n;
+      int64_t k64=k;
+      int64_t lda64=lda;
+      int64_t ldb64=ldb;
+      int64_t ldc64=ldc;
+
+      oneapi::mkl::transpose iOpA;
+      oneapi::mkl::transpose iOpB;
+      
+      if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
+      if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
+      if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
+      if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
+      if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
+      if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
+
+      oneapi::mkl::blas::column_major::gemm(*gridblasHandle,
+					    iOpA,
+					    iOpB,
+					    m64,n64,k64,
+					    (ComplexF *) &alpha_p[0],
+					    (const ComplexF *)Amk, (int64_t )lda64,
+					    (const ComplexF *)Bkn, (int64_t )ldb64,
+					    (ComplexF *) &beta_p[0],
+					    (ComplexF *)Cmn, (int64_t )ldc64);
+      synchronise();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
+    // Need a default/reference implementation; use Eigen
+      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+	Eigen::Map<Eigen::MatrixXcf> eAmk(Amk,m,k);
+	Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn,k,n);
+	Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+	Eigen::Map<Eigen::MatrixXcf> eAmk(Amk,k,m);
+	Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn,k,n);
+	Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn ;
+      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
+	Eigen::Map<Eigen::MatrixXcf> eAmk(Amk,m,k);
+	Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn,n,k);
+	Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint() ;
+      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
+	Eigen::Map<Eigen::MatrixXcf> eAmk(Amk,k,m);
+	Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn,n,k);
+	Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn,m,n);
+	eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint() ;
+      } else { 
+	assert(0);
+      }
+#endif
+     RealD t1=usecond();
+     RealD flops = 8.0*m*n*k;
+     RealD bytes = 1.0*sizeof(ComplexF)*(m*k+k*n+m*n);
+  }
+
   
+  /////////////////////////////////////////////////////////////
   void gemmBatched(int m,int n, int k,
 		   ComplexD alpha,
 		   deviceVector<ComplexD*> &Amk,  // pointer list to matrices
@@ -241,36 +583,6 @@ public:
 		beta,
 		Cmn);
   }
-  void gemmBatched(int m,int n, int k,
-		   RealD alpha,
-		   deviceVector<RealD*> &Amk,  // pointer list to matrices
-		   deviceVector<RealD*> &Bkn,
-		   RealD beta,
-		   deviceVector<RealD*> &Cmn)
-  {
-    gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
-		m,n,k,
-		alpha,
-		Amk,
-		Bkn,
-		beta,
-		Cmn);
-  }
-  void gemmBatched(int m,int n, int k,
-		   RealF alpha,
-		   deviceVector<RealF*> &Amk,  // pointer list to matrices
-		   deviceVector<RealF*> &Bkn,
-		   RealF beta,
-		   deviceVector<RealF*> &Cmn)
-  {
-    gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
-		m,n,k,
-		alpha,
-		Amk,
-		Bkn,
-		beta,
-		Cmn);
-  }
 
   void gemmBatched(GridBLASOperation_t OpA,
 		   GridBLASOperation_t OpB,
@@ -623,301 +935,6 @@ public:
      RealD flops = 8.0*m*n*k*batchCount;
      RealD bytes = 1.0*sizeof(ComplexF)*(m*k+k*n+m*n)*batchCount;
   }
-  
-  ///////////////////////////////////////////////////////////////////////////
-  // Single precision real GEMM
-  ///////////////////////////////////////////////////////////////////////////
-
-  void gemmBatched(GridBLASOperation_t OpA,
-		   GridBLASOperation_t OpB,
-		   int m,int n, int k,
-		   RealF alpha,
-		   deviceVector<RealF*> &Amk,  // pointer list to matrices
-		   deviceVector<RealF*> &Bkn,
-		   RealF beta,
-		   deviceVector<RealF*> &Cmn)
-  {
-    RealD t2=usecond();
-    int32_t batchCount = Amk.size();
-
-    assert(OpA!=GridBLAS_OP_C); // Real case no conjugate
-    assert(OpB!=GridBLAS_OP_C);
-
-    int lda = m; // m x k column major
-    int ldb = k; // k x n column major
-    int ldc = m; // m x b column major
-    if(OpA!=GridBLAS_OP_N)
-      lda = k;
-    if(OpB!=GridBLAS_OP_N)
-      ldb = n;
-    static deviceVector<RealF> alpha_p(1);
-    static deviceVector<RealF> beta_p(1);
-    // can prestore the 1 and the zero on device
-    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(RealF));
-    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(RealF));
-    RealD t0=usecond();
-
-    assert(Bkn.size()==batchCount);
-    assert(Cmn.size()==batchCount);
-#ifdef GRID_HIP
-    hipblasOperation_t hOpA;
-    hipblasOperation_t hOpB;
-    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
-    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
-    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
-    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
-    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
-    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
-    auto err = hipblasSgemmBatched(gridblasHandle,
-				   hOpA,
-				   hOpB,
-				   m,n,k,
-				   (float *) &alpha_p[0],
-				   (float **)&Amk[0], lda,
-				   (float **)&Bkn[0], ldb,
-				   (float *) &beta_p[0],
-				   (float **)&Cmn[0], ldc,
-				   batchCount);
-    assert(err==HIPBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_CUDA
-    cublasOperation_t hOpA;
-    cublasOperation_t hOpB;
-    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
-    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
-    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
-    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
-    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
-    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
-    auto err = cublasSgemmBatched(gridblasHandle,
-				  hOpA,
-				  hOpB,
-				  m,n,k,
-				  (float *) &alpha_p[0],
-				  (float **)&Amk[0], lda,
-				  (float **)&Bkn[0], ldb,
-				  (float *) &beta_p[0],
-				  (float **)&Cmn[0], ldc,
-				  batchCount);
-    assert(err==CUBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_SYCL
-      int64_t m64=m;
-      int64_t n64=n;
-      int64_t k64=k;
-      int64_t lda64=lda;
-      int64_t ldb64=ldb;
-      int64_t ldc64=ldc;
-      int64_t batchCount64=batchCount;
-
-      oneapi::mkl::transpose iOpA;
-      oneapi::mkl::transpose iOpB;
-      
-      if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
-      if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
-      if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
-      if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
-      if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
-      if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
-
-      oneapi::mkl::blas::column_major::gemm_batch(*gridblasHandle,
-						  &iOpA,
-						  &iOpB,
-						  &m64,&n64,&k64,
-						  (float *) &alpha_p[0],
-						  (const float **)&Amk[0], (const int64_t *)&lda64,
-						  (const float **)&Bkn[0], (const int64_t *)&ldb64,
-						  (float *) &beta_p[0],
-						  (float **)&Cmn[0], (const int64_t *)&ldc64,
-						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
-      synchronise();
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXf> eBkn(Bkn[p],k,n);
-	  Eigen::Map<Eigen::MatrixXf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXf> eBkn(Bkn[p],k,n);
-	  Eigen::Map<Eigen::MatrixXf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_T) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXf> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn.transpose() ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_T) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXf> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn.transpose() ;
-	  } );
-      } else { 
-	assert(0);
-      }
-#endif
-     RealD t1=usecond();
-     RealD flops = 2.0*m*n*k*batchCount;
-     RealD bytes = 1.0*sizeof(RealF)*(m*k+k*n+m*n)*batchCount;
-  }
-  
-  
-  ///////////////////////////////////////////////////////////////////////////
-  // Double precision real GEMM
-  ///////////////////////////////////////////////////////////////////////////
-  void gemmBatched(GridBLASOperation_t OpA,
-		   GridBLASOperation_t OpB,
-		   int m,int n, int k,
-		   RealD alpha,
-		   deviceVector<RealD*> &Amk,  // pointer list to matrices
-		   deviceVector<RealD*> &Bkn,
-		   RealD beta,
-		   deviceVector<RealD*> &Cmn)
-  {
-    RealD t2=usecond();
-    int32_t batchCount = Amk.size();
-
-    assert(OpA!=GridBLAS_OP_C); // Real case no conjugate
-    assert(OpB!=GridBLAS_OP_C);
-
-    int lda = m; // m x k column major
-    int ldb = k; // k x n column major
-    int ldc = m; // m x b column major
-    if(OpA!=GridBLAS_OP_N)
-      lda = k;
-    if(OpB!=GridBLAS_OP_N)
-      ldb = n;
-    
-    static deviceVector<RealD> alpha_p(1);
-    static deviceVector<RealD> beta_p(1);
-    // can prestore the 1 and the zero on device
-    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(RealD));
-    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(RealD));
-    RealD t0=usecond();
-
-    assert(Bkn.size()==batchCount);
-    assert(Cmn.size()==batchCount);
-#ifdef GRID_HIP
-    hipblasOperation_t hOpA;
-    hipblasOperation_t hOpB;
-    if ( OpA == GridBLAS_OP_N ) hOpA = HIPBLAS_OP_N;
-    if ( OpA == GridBLAS_OP_T ) hOpA = HIPBLAS_OP_T;
-    if ( OpA == GridBLAS_OP_C ) hOpA = HIPBLAS_OP_C;
-    if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
-    if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
-    if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
-    auto err = hipblasDgemmBatched(gridblasHandle,
-				   HIPBLAS_OP_N,
-				   HIPBLAS_OP_N,
-				   m,n,k,
-				   (double *) &alpha_p[0],
-				   (double **)&Amk[0], lda,
-				   (double **)&Bkn[0], ldb,
-				   (double *) &beta_p[0],
-				   (double **)&Cmn[0], ldc,
-				   batchCount);
-    assert(err==HIPBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_CUDA
-    cublasOperation_t hOpA;
-    cublasOperation_t hOpB;
-    if ( OpA == GridBLAS_OP_N ) hOpA = CUBLAS_OP_N;
-    if ( OpA == GridBLAS_OP_T ) hOpA = CUBLAS_OP_T;
-    if ( OpA == GridBLAS_OP_C ) hOpA = CUBLAS_OP_C;
-    if ( OpB == GridBLAS_OP_N ) hOpB = CUBLAS_OP_N;
-    if ( OpB == GridBLAS_OP_T ) hOpB = CUBLAS_OP_T;
-    if ( OpB == GridBLAS_OP_C ) hOpB = CUBLAS_OP_C;
-    auto err = cublasDgemmBatched(gridblasHandle,
-				  hOpA,
-				  hOpB,
-				  m,n,k,
-				  (double *) &alpha_p[0],
-				  (double **)&Amk[0], lda,
-				  (double **)&Bkn[0], ldb,
-				  (double *) &beta_p[0],
-				  (double **)&Cmn[0], ldc,
-				  batchCount);
-    assert(err==CUBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_SYCL
-      int64_t m64=m;
-      int64_t n64=n;
-      int64_t k64=k;
-      int64_t lda64=lda;
-      int64_t ldb64=ldb;
-      int64_t ldc64=ldc;
-      int64_t batchCount64=batchCount;
-
-      oneapi::mkl::transpose iOpA;
-      oneapi::mkl::transpose iOpB;
-      
-      if ( OpA == GridBLAS_OP_N ) iOpA = oneapi::mkl::transpose::N;
-      if ( OpA == GridBLAS_OP_T ) iOpA = oneapi::mkl::transpose::T;
-      if ( OpA == GridBLAS_OP_C ) iOpA = oneapi::mkl::transpose::C;
-      if ( OpB == GridBLAS_OP_N ) iOpB = oneapi::mkl::transpose::N;
-      if ( OpB == GridBLAS_OP_T ) iOpB = oneapi::mkl::transpose::T;
-      if ( OpB == GridBLAS_OP_C ) iOpB = oneapi::mkl::transpose::C;
-
-      oneapi::mkl::blas::column_major::gemm_batch(*gridblasHandle,
-						  &iOpA,
-						  &iOpB,
-						  &m64,&n64,&k64,
-						  (double *) &alpha_p[0],
-						  (const double **)&Amk[0], (const int64_t *)&lda64,
-						  (const double **)&Bkn[0], (const int64_t *)&ldb64,
-						  (double *) &beta_p[0],
-						  (double **)&Cmn[0], (const int64_t *)&ldc64,
-						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
-      synchronise();
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXd> eBkn(Bkn[p],k,n);
-	  Eigen::Map<Eigen::MatrixXd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXd> eBkn(Bkn[p],k,n);
-	  Eigen::Map<Eigen::MatrixXd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_T) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXd> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn.transpose() ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_T) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXd> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.transpose() * eBkn.transpose() ;
-	  });
-      } else { 
-	assert(0);
-      }
-#endif
-     RealD t1=usecond();
-     RealD flops = 2.0*m*n*k*batchCount;
-     RealD bytes = 1.0*sizeof(RealD)*(m*k+k*n+m*n)*batchCount;
-  }
 
   template<class CComplex>
   double benchmark(int M, int N, int K, int BATCH)
@@ -967,6 +984,47 @@ public:
     return flops; // Returns gigaflops
   }
 
+  template<class CComplex>
+  double benchmark(int M, int N, int K)
+  {
+    int32_t N_A = M*K;
+    int32_t N_B = K*N;
+    int32_t N_C = M*N;
+    deviceVector<CComplex> A(N_A); acceleratorMemSet(&A[0],0,N_A*sizeof(CComplex));
+    deviceVector<CComplex> B(N_B); acceleratorMemSet(&B[0],0,N_B*sizeof(CComplex));
+    deviceVector<CComplex> C(N_C); acceleratorMemSet(&C[0],0,N_C*sizeof(CComplex));
+    CComplex alpha(1.0);
+    CComplex beta (1.0);
+    RealD flops = 8.0*M*N*K;
+    int ncall=10;
+
+    gemm(GridBLAS_OP_C,GridBLAS_OP_N,
+	 M,N,K,
+	 alpha,
+	 &A[0], // m x k 
+	 &B[0], // k x n
+	 beta, 
+	 &C[0]);
+    synchronise();
+
+    RealD t0 = usecond();
+    for(int i=0;i<ncall;i++){
+      gemm(GridBLAS_OP_N,GridBLAS_OP_N,
+	   M,N,K,
+	   alpha,
+	   &A[0], // m x k 
+	   &B[0], // k x n
+	   beta, 
+	   &C[0]);
+      synchronise();
+    }
+    RealD t1 = usecond();
+    RealD bytes = 1.0*sizeof(CComplex)*(M*N*2+N*K+M*K);
+    flops = 8.0*M*N*K*ncall;
+    flops = flops/(t1-t0)/1.e3;
+    return flops; // Returns gigaflops
+  }
+
 };
 
 
@@ -1035,6 +1093,21 @@ static void BLAS(void)
       std::cout<< M<<"\t\t"<<N<<"\t\t"<<K<<"\t\t"<<BATCH<<"\t\t"<<p<<std::endl;
     }}
   fprintf(FP,"\n\n\n");
+
+  std::cout << "----------------------------------------------------------"<<std::endl;
+  std::cout << "  M  "<<"\t\t"<<"N"<<"\t\t\t"<<"K"<<"\t\t"<<"Gflop/s / rank (inner product matrix)"<<std::endl;
+  std::cout << "----------------------------------------------------------"<<std::endl;
+  {
+    int M=12;
+    int N=12;
+    std::vector<int> ks({4*1024*1024, 2*1024*1024, 1024*1024, 256*1024, 1024 });
+    for( int kk=0;kk<ks.size();kk++ ) {
+      int K = ks[kk];
+      double p=blas.benchmark<CComplex>(M,N,K);
+      fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, 1, p);
+      std::cout<< M<<"\t\t"<<N<<"\t\t"<<K<<"\t\t"<<1<<"\t\t"<<p<<std::endl;
+    }
+  }
   std::cout << "=================================================================================="<<std::endl;
 };
 

BLAS_benchmark/compile-command

@@ -1,2 +1,2 @@
 
-mpicxx -qmkl=parallel -fsycl BatchBlasBench.cc -o BatchBlasBench
+mpicxx -qmkl=parallel -fsycl BatchBlasBench.cc -o BatchBlasBench -DGRID_SYCL

BLAS_benchmark/compile-command-frontier

@@ -0,0 +1,5 @@
+CXX=hipcc
+MPICXX=mpicxx 
+CXXFLAGS="-fPIC -I{$ROCM_PATH}/include/ -I${MPICH_DIR}/include -L/lib64 -I/opt/cray/pe/mpich/8.1.28/ofi/gnu/12.3/include -DGRID_HIP"
+LDFLAGS="-L/lib64 -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64 -lhipblas -lrocblas -lmpi_gnu_123"
+hipcc $CXXFLAGS $LDFLAGS BatchBlasBench.cc -o BatchBlasBench

BLAS_benchmark/compile-command-sunspot

@@ -0,0 +1,2 @@
+
+mpicxx -qmkl=parallel -fsycl BatchBlasBench.cc -o BatchBlasBench -DGRID_SYCL

Grid/algorithms/Algorithms.h

@@ -50,6 +50,7 @@ NAMESPACE_CHECK(approx);
 #include <Grid/algorithms/deflation/Deflation.h>
 #include <Grid/algorithms/deflation/MultiRHSBlockProject.h>
 #include <Grid/algorithms/deflation/MultiRHSDeflation.h>
+#include <Grid/algorithms/deflation/MultiRHSBlockCGLinalg.h>
 NAMESPACE_CHECK(deflation);
 #include <Grid/algorithms/iterative/ConjugateGradient.h>
 NAMESPACE_CHECK(ConjGrad);

Grid/algorithms/LinearOperator.h

@@ -103,6 +103,38 @@ public:
     _Mat.MdagM(in,out);
   }
 };
+template<class Matrix,class Field>
+class MMdagLinearOperator : public LinearOperatorBase<Field> {
+  Matrix &_Mat;
+public:
+  MMdagLinearOperator(Matrix &Mat): _Mat(Mat){};
+
+  // Support for coarsening to a multigrid
+  void OpDiag (const Field &in, Field &out) {
+    _Mat.Mdiag(in,out);
+  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {
+    _Mat.Mdir(in,out,dir,disp);
+  }
+  void OpDirAll  (const Field &in, std::vector<Field> &out){
+    _Mat.MdirAll(in,out);
+  };
+  void Op     (const Field &in, Field &out){
+    _Mat.M(in,out);
+  }
+  void AdjOp     (const Field &in, Field &out){
+    _Mat.Mdag(in,out);
+  }
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
+    _Mat.MMdag(in,out);
+    ComplexD dot = innerProduct(in,out);
+    n1=real(dot);
+    n2=norm2(out);
+  }
+  void HermOp(const Field &in, Field &out){
+    _Mat.MMdag(in,out);
+  }
+};
 
 ////////////////////////////////////////////////////////////////////
 // Construct herm op and shift it for mgrid smoother

Grid/algorithms/SparseMatrix.h

@@ -45,6 +45,11 @@ public:
     M(in,tmp);
     Mdag(tmp,out);
   }
+  virtual void  MMdag(const Field &in, Field &out) {
+    Field tmp (in.Grid());
+    Mdag(in,tmp);
+    M(tmp,out);
+  }
   virtual  void Mdiag    (const Field &in, Field &out)=0;
   virtual  void Mdir     (const Field &in, Field &out,int dir, int disp)=0;
   virtual  void MdirAll  (const Field &in, std::vector<Field> &out)=0;

Grid/algorithms/approx/Chebyshev.h

@@ -59,7 +59,7 @@ public:
     RealD diff = hi-lo;
     RealD delta = diff*1.0e-9;
     for (RealD x=lo; x<hi; x+=delta) {
-      delta*=1.1;
+      delta*=1.02;
       RealD f = approx(x);
       out<< x<<" "<<f<<std::endl;
     }
@@ -131,6 +131,26 @@ public:
       Coeffs[j] = s * 2.0/order;
     }
   };
+  template<class functor>
+  void Init(RealD _lo,RealD _hi,int _order, functor & func)
+  {
+    lo=_lo;
+    hi=_hi;
+    order=_order;
+      
+    if(order < 2) exit(-1);
+    Coeffs.resize(order);
+    for(int j=0;j<order;j++){
+      RealD s=0;
+      for(int k=0;k<order;k++){
+	RealD y=std::cos(M_PI*(k+0.5)/order);
+	RealD x=0.5*(y*(hi-lo)+(hi+lo));
+	RealD f=func(x);
+	s=s+f*std::cos( j*M_PI*(k+0.5)/order );
+      }
+      Coeffs[j] = s * 2.0/order;
+    }
+  };
 
     
   void JacksonSmooth(void){

Grid/algorithms/deflation/MultiRHSBlockCGLinalg.h

@@ -0,0 +1,376 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: MultiRHSBlockCGLinalg.h
+
+    Copyright (C) 2024
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+
+/* Need helper object for BLAS accelerated mrhs blockCG */
+template<class Field>
+class MultiRHSBlockCGLinalg
+{
+public:
+
+  typedef typename Field::scalar_type   scalar;
+  typedef typename Field::scalar_object scalar_object;
+  typedef typename Field::vector_object vector_object;
+
+  deviceVector<scalar> BLAS_X;      // nrhs x vol -- the sources
+  deviceVector<scalar> BLAS_Y;      // nrhs x vol -- the result
+  deviceVector<scalar> BLAS_C;      // nrhs x nrhs -- the coefficients 
+  deviceVector<scalar> BLAS_Cred;   // nrhs x nrhs x oSites -- reduction buffer
+  deviceVector<scalar *> Xdip;
+  deviceVector<scalar *> Ydip;
+  deviceVector<scalar *> Cdip;
+  
+  MultiRHSBlockCGLinalg() {};
+  ~MultiRHSBlockCGLinalg(){ Deallocate(); };
+  
+  void Deallocate(void)
+  {
+    Xdip.resize(0);
+    Ydip.resize(0);
+    Cdip.resize(0);
+    BLAS_Cred.resize(0);
+    BLAS_C.resize(0);
+    BLAS_X.resize(0);
+    BLAS_Y.resize(0);
+  }
+  void MaddMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X,const std::vector<Field> &Y,RealD scale=1.0)
+  {
+    std::vector<Field> Y_copy(AP.size(),AP[0].Grid());
+    for(int r=0;r<AP.size();r++){
+      Y_copy[r] = Y[r];
+    }
+    MulMatrix(AP,m,X);
+    for(int r=0;r<AP.size();r++){
+      AP[r] = scale*AP[r]+Y_copy[r];
+    }
+  }
+  void MulMatrix(std::vector<Field> &Y, Eigen::MatrixXcd &m , const std::vector<Field> &X)
+  {
+    typedef typename Field::scalar_type scomplex;
+    GridBase *grid;
+    uint64_t vol;
+    uint64_t words;
+
+    int nrhs = Y.size();
+    grid  = X[0].Grid();
+    vol   = grid->lSites();
+    words = sizeof(scalar_object)/sizeof(scalar);
+    int64_t vw = vol * words;
+
+    RealD t0 = usecond();
+    BLAS_X.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_Y.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_C.resize(nrhs * nrhs);// cost free if size doesn't change
+    RealD t1 = usecond();
+
+    /////////////////////////////////////////////
+    // Copy in the multi-rhs sources
+    /////////////////////////////////////////////
+    for(int r=0;r<nrhs;r++){
+      int64_t offset = r*vw;
+      autoView(x_v,X[r],AcceleratorRead);
+      acceleratorCopyDeviceToDevice(&x_v[0],&BLAS_X[offset],sizeof(scalar_object)*vol);
+    }
+
+    // Assumes Eigen storage contiguous
+    acceleratorCopyToDevice(&m(0,0),&BLAS_C[0],BLAS_C.size()*sizeof(scalar));
+    
+  /*
+   * in Fortran column major notation (cuBlas order)
+   *
+   * Xxr = [X1(x)][..][Xn(x)]
+   * Yxr = [Y1(x)][..][Ym(x)]
+   * Y = X . C
+   */
+    deviceVector<scalar *> Xd(1);
+    deviceVector<scalar *> Yd(1);
+    deviceVector<scalar *> Cd(1);
+
+    scalar * Xh = & BLAS_X[0];
+    scalar * Yh = & BLAS_Y[0];
+    scalar * Ch = & BLAS_C[0];
+
+    acceleratorPut(Xd[0],Xh);
+    acceleratorPut(Yd[0],Yh);
+    acceleratorPut(Cd[0],Ch);
+
+    RealD t2 = usecond();
+    GridBLAS BLAS;
+    /////////////////////////////////////////
+    // Y = X*C (transpose?)
+    /////////////////////////////////////////
+    BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N, 
+    		     vw,nrhs,nrhs,
+		     scalar(1.0),
+		     Xd,
+		     Cd,
+		     scalar(0.0),  // wipe out Y
+		     Yd);
+    BLAS.synchronise();
+    RealD t3 = usecond();
+
+    // Copy back Y = m X 
+    for(int r=0;r<nrhs;r++){
+      int64_t offset = r*vw;
+      autoView(y_v,Y[r],AcceleratorWrite);
+      acceleratorCopyDeviceToDevice(&BLAS_Y[offset],&y_v[0],sizeof(scalar_object)*vol);
+    }    
+    RealD t4 = usecond();
+    std::cout << "MulMatrix alloc    took "<< t1-t0<<" us"<<std::endl;
+    std::cout << "MulMatrix preamble took "<< t2-t1<<" us"<<std::endl;
+    std::cout << "MulMatrix blas     took "<< t3-t2<<" us"<<std::endl;
+    std::cout << "MulMatrix copy     took "<< t4-t3<<" us"<<std::endl;
+    std::cout << "MulMatrix total "<< t4-t0<<" us"<<std::endl;
+  }
+  
+  void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector<Field> &X, const std::vector<Field> &Y)
+  {
+#if 0    
+    int nrhs;
+    GridBase *grid;
+    uint64_t vol;
+    uint64_t words;
+
+    nrhs = X.size();
+    assert(X.size()==Y.size());
+    conformable(X[0],Y[0]);
+
+    grid  = X[0].Grid();
+    vol   = grid->lSites();
+    words = sizeof(scalar_object)/sizeof(scalar);
+    int64_t vw = vol * words;
+
+    RealD t0 = usecond();
+    BLAS_X.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_Y.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_C.resize(nrhs * nrhs);// cost free if size doesn't change
+    RealD t1 = usecond();
+
+    /////////////////////////////////////////////
+    // Copy in the multi-rhs sources
+    /////////////////////////////////////////////
+    for(int r=0;r<nrhs;r++){
+      int64_t offset = r*vw;
+      autoView(x_v,X[r],AcceleratorRead);
+      acceleratorCopyDeviceToDevice(&x_v[0],&BLAS_X[offset],sizeof(scalar_object)*vol);
+      autoView(y_v,Y[r],AcceleratorRead);
+      acceleratorCopyDeviceToDevice(&y_v[0],&BLAS_Y[offset],sizeof(scalar_object)*vol);
+    }
+    RealD t2 = usecond();
+
+  /*
+   * in Fortran column major notation (cuBlas order)
+   *
+   * Xxr = [X1(x)][..][Xn(x)]
+   *
+   * Yxr = [Y1(x)][..][Ym(x)]
+   *
+   * C_rs = X^dag Y
+   */
+    deviceVector<scalar *> Xd(1);
+    deviceVector<scalar *> Yd(1);
+    deviceVector<scalar *> Cd(1);
+
+    scalar * Xh = & BLAS_X[0];
+    scalar * Yh = & BLAS_Y[0];
+    scalar * Ch = & BLAS_C[0];
+
+    acceleratorPut(Xd[0],Xh);
+    acceleratorPut(Yd[0],Yh);
+    acceleratorPut(Cd[0],Ch);
+
+    GridBLAS BLAS;
+
+    RealD t3 = usecond();
+    /////////////////////////////////////////
+    // C_rs = X^dag Y
+    /////////////////////////////////////////
+    BLAS.gemmBatched(GridBLAS_OP_C,GridBLAS_OP_N, 
+    		     nrhs,nrhs,vw,
+		     ComplexD(1.0),
+		     Xd,
+		     Yd,
+		     ComplexD(0.0),  // wipe out C
+		     Cd);
+    BLAS.synchronise();
+    RealD t4 = usecond();
+
+    std::vector<scalar> HOST_C(BLAS_C.size());      // nrhs . nrhs -- the coefficients 
+    acceleratorCopyFromDevice(&BLAS_C[0],&HOST_C[0],BLAS_C.size()*sizeof(scalar));
+    grid->GlobalSumVector(&HOST_C[0],nrhs*nrhs);
+
+    RealD t5 = usecond();
+    for(int rr=0;rr<nrhs;rr++){
+      for(int r=0;r<nrhs;r++){
+	int off = r+nrhs*rr;
+	m(r,rr)=HOST_C[off];
+      }
+    }
+    RealD t6 = usecond();
+    uint64_t M=nrhs;
+    uint64_t N=nrhs;
+    uint64_t K=vw;
+    RealD bytes = 1.0*sizeof(ComplexD)*(M*N*2+N*K+M*K);
+    RealD flops = 8.0*M*N*K;
+    flops = flops/(t4-t3)/1.e3;
+    bytes = bytes/(t4-t3)/1.e3;
+    std::cout << "InnerProductMatrix m,n,k "<< M<<","<<N<<","<<K<<std::endl;
+    std::cout << "InnerProductMatrix alloc t1 "<< t1-t0<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix cp    t2 "<< t2-t1<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix setup t3 "<< t3-t2<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix blas t4 "<< t4-t3<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix blas    "<< flops<<" GF/s"<<std::endl;
+    std::cout << "InnerProductMatrix blas    "<< bytes<<" GB/s"<<std::endl;
+    std::cout << "InnerProductMatrix gsum t5 "<< t5-t4<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix cp   t6 "<< t6-t5<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix took "<< t6-t0<<" us"<<std::endl;
+#else
+    int nrhs;
+    GridBase *grid;
+    uint64_t vol;
+    uint64_t words;
+
+    nrhs = X.size();
+    assert(X.size()==Y.size());
+    conformable(X[0],Y[0]);
+
+    grid  = X[0].Grid();
+    int rd0 =  grid->_rdimensions[0] * grid->_rdimensions[1];
+    vol   = grid->oSites()/rd0;
+    words = rd0*sizeof(vector_object)/sizeof(scalar);
+    int64_t vw = vol * words;
+    assert(vw == grid->lSites()*sizeof(scalar_object)/sizeof(scalar));
+
+    RealD t0 = usecond();
+    BLAS_X.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_Y.resize(nrhs * vw); // cost free if size doesn't change
+    BLAS_Cred.resize(nrhs * nrhs * vol);// cost free if size doesn't change
+    RealD t1 = usecond();
+
+    /////////////////////////////////////////////
+    // Copy in the multi-rhs sources -- layout batched BLAS ready
+    /////////////////////////////////////////////
+    for(int r=0;r<nrhs;r++){
+      autoView(x_v,X[r],AcceleratorRead);
+      autoView(y_v,Y[r],AcceleratorRead);
+      scalar *from_x=(scalar *)&x_v[0];
+      scalar *from_y=(scalar *)&y_v[0];
+      scalar *BX = &BLAS_X[0];
+      scalar *BY = &BLAS_Y[0];
+      accelerator_for(ssw,vw,1,{
+	  uint64_t ss=ssw/words;
+	  uint64_t  w=ssw%words;
+	  uint64_t offset = w+r*words+ss*nrhs*words; // [ss][rhs][words]
+	  BX[offset] = from_x[ssw];
+	  BY[offset] = from_y[ssw];
+	});
+    }
+    RealD t2 = usecond();
+
+  /*
+   * in Fortran column major notation (cuBlas order)
+   *
+   * Xxr = [X1(x)][..][Xn(x)]
+   *
+   * Yxr = [Y1(x)][..][Ym(x)]
+   *
+   * C_rs = X^dag Y
+   */
+    Xdip.resize(vol);
+    Ydip.resize(vol);
+    Cdip.resize(vol);
+    std::vector<scalar *> Xh(vol);
+    std::vector<scalar *> Yh(vol);
+    std::vector<scalar *> Ch(vol);
+    for(uint64_t ss=0;ss<vol;ss++){
+
+      Xh[ss] = & BLAS_X[ss*nrhs*words];
+      Yh[ss] = & BLAS_Y[ss*nrhs*words];
+      Ch[ss] = & BLAS_Cred[ss*nrhs*nrhs];
+
+    }
+    acceleratorCopyToDevice(&Xh[0],&Xdip[0],vol*sizeof(scalar *));
+    acceleratorCopyToDevice(&Yh[0],&Ydip[0],vol*sizeof(scalar *));
+    acceleratorCopyToDevice(&Ch[0],&Cdip[0],vol*sizeof(scalar *));
+    
+    GridBLAS BLAS;
+
+    RealD t3 = usecond();
+    /////////////////////////////////////////
+    // C_rs = X^dag Y
+    /////////////////////////////////////////
+    BLAS.gemmBatched(GridBLAS_OP_C,GridBLAS_OP_N, 
+    		     nrhs,nrhs,words,
+		     ComplexD(1.0),
+		     Xdip,
+		     Ydip,
+		     ComplexD(0.0),  // wipe out C
+		     Cdip);
+    BLAS.synchronise();
+    RealD t4 = usecond();
+
+    std::vector<scalar> HOST_C(BLAS_Cred.size());      // nrhs . nrhs -- the coefficients 
+    acceleratorCopyFromDevice(&BLAS_Cred[0],&HOST_C[0],BLAS_Cred.size()*sizeof(scalar));
+
+    RealD t5 = usecond();
+    m = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    for(int ss=0;ss<vol;ss++){
+      Eigen::Map<Eigen::MatrixXcd> eC((std::complex<double> *)&HOST_C[ss*nrhs*nrhs],nrhs,nrhs);
+      m = m + eC;
+    }
+    RealD t6l = usecond();
+    grid->GlobalSumVector((scalar *) &m(0,0),nrhs*nrhs);
+    RealD t6 = usecond();
+    uint64_t M=nrhs;
+    uint64_t N=nrhs;
+    uint64_t K=vw;
+    RealD xybytes = grid->lSites()*sizeof(scalar_object);
+    RealD bytes = 1.0*sizeof(ComplexD)*(M*N*2+N*K+M*K);
+    RealD flops = 8.0*M*N*K;
+    flops = flops/(t4-t3)/1.e3;
+    bytes = bytes/(t4-t3)/1.e3;
+    xybytes = 4*xybytes/(t2-t1)/1.e3;
+    std::cout << "InnerProductMatrix m,n,k "<< M<<","<<N<<","<<K<<std::endl;
+    std::cout << "InnerProductMatrix alloc t1 "<< t1-t0<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix cp    t2 "<< t2-t1<<" us "<<xybytes<<" GB/s"<<std::endl;
+    std::cout << "InnerProductMatrix setup t3 "<< t3-t2<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix blas t4 "<< t4-t3<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix blas    "<< flops<<" GF/s"<<std::endl;
+    std::cout << "InnerProductMatrix blas    "<< bytes<<" GB/s"<<std::endl;
+    std::cout << "InnerProductMatrix cp     t5 "<< t5-t4<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix lsum   t6l "<< t6l-t5<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix gsum   t6 "<< t6-t6l<<" us"<<std::endl;
+    std::cout << "InnerProductMatrix took "<< t6-t0<<" us"<<std::endl;
+#endif
+  }
+};
+
+NAMESPACE_END(Grid);

Grid/algorithms/deflation/MultiRHSBlockProject.h

@@ -447,10 +447,10 @@ public:
     /////////////////////////////////////////
     BLAS.gemmBatched(GridBLAS_OP_C,GridBLAS_OP_N, 
     		     nbasis,nrhs,vw,
-		     ComplexD(1.0),
+		     scalar(1.0),
 		     Vd,
 		     Fd,
-		     ComplexD(0.0),  // wipe out C
+		     scalar(0.0),  // wipe out C
 		     Cd);
     BLAS.synchronise();
     //    std::cout << "BlockProject done"<<std::endl;
@@ -497,10 +497,10 @@ public:
     int64_t vw = block_vol * words;
     BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N, 
     		     vw,nrhs,nbasis,
-		     ComplexD(1.0),
+		     scalar(1.0),
 		     Vd,
 		     Cd,
-		     ComplexD(0.0),  // wipe out C
+		     scalar(0.0),  // wipe out C
 		     Fd);
     BLAS.synchronise();
     //    std::cout << " blas call done"<<std::endl;

Grid/algorithms/deflation/MultiRHSDeflation.h

@@ -182,10 +182,10 @@ public:
     /////////////////////////////////////////
     BLAS.gemmBatched(GridBLAS_OP_C,GridBLAS_OP_N, 
     		     nev,nrhs,vw,
-		     ComplexD(1.0),
+		     scalar(1.0),
 		     Ed,
 		     Rd,
-		     ComplexD(0.0),  // wipe out C
+		     scalar(0.0),  // wipe out C
 		     Cd);
     BLAS.synchronise();
 
@@ -210,10 +210,10 @@ public:
     /////////////////////////////////////////
     BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N, 
 		     vw,nrhs,nev,
-		     ComplexD(1.0),
+		     scalar(1.0),
 		     Ed, // x . nev
 		     Cd, // nev . nrhs
-		     ComplexD(0.0),
+		     scalar(0.0),
 		     Gd);
     BLAS.synchronise();
 

Grid/algorithms/iterative/AdefMrhs.h

@@ -53,6 +53,7 @@ class TwoLevelCGmrhs
   // Fine operator, Smoother, CoarseSolver
   LinearOperatorBase<Field>   &_FineLinop;
   LinearFunction<Field>   &_Smoother;
+  MultiRHSBlockCGLinalg<Field> _BlockCGLinalg;
 
   GridStopWatch ProjectTimer;
   GridStopWatch PromoteTimer;
@@ -79,6 +80,301 @@ class TwoLevelCGmrhs
   
   // Vector case
   virtual void operator() (std::vector<Field> &src, std::vector<Field> &x)
+  {
+    SolveSingleSystem(src,x);
+    //    SolvePrecBlockCG(src,x);
+  }
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Thin QR factorisation (google it)
+////////////////////////////////////////////////////////////////////////////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////////////////////////////
+  //Dimensions
+  // R_{ferm x Nblock} =  Q_{ferm x Nblock} x  C_{Nblock x Nblock} -> ferm x Nblock
+  //
+  // Rdag R = m_rr = Herm = L L^dag        <-- Cholesky decomposition (LLT routine in Eigen)
+  //
+  //   Q  C = R => Q = R C^{-1}
+  //
+  // Want  Ident = Q^dag Q = C^{-dag} R^dag R C^{-1} = C^{-dag} L L^dag C^{-1} = 1_{Nblock x Nblock} 
+  //
+  // Set C = L^{dag}, and then Q^dag Q = ident 
+  //
+  // Checks:
+  // Cdag C = Rdag R ; passes.
+  // QdagQ  = 1      ; passes
+  ////////////////////////////////////////////////////////////////////////////////////////////////////
+  void ThinQRfact (Eigen::MatrixXcd &m_zz,
+		   Eigen::MatrixXcd &C,
+		   Eigen::MatrixXcd &Cinv,
+		   std::vector<Field> &  Q,
+		   std::vector<Field> & MQ,
+		   const std::vector<Field> & Z,
+		   const std::vector<Field> & MZ)
+  {
+    RealD t0=usecond();
+    _BlockCGLinalg.InnerProductMatrix(m_zz,MZ,Z);
+    RealD t1=usecond();
+
+    m_zz = 0.5*(m_zz+m_zz.adjoint());
+    
+    Eigen::MatrixXcd L    = m_zz.llt().matrixL(); 
+    
+    C    = L.adjoint();
+    Cinv = C.inverse();
+    
+    RealD t3=usecond();
+    _BlockCGLinalg.MulMatrix( Q,Cinv,Z);
+    _BlockCGLinalg.MulMatrix(MQ,Cinv,MZ);
+    RealD t4=usecond();
+    std::cout << " ThinQRfact IP    :"<< t1-t0<<" us"<<std::endl;
+    std::cout << " ThinQRfact Eigen :"<< t3-t1<<" us"<<std::endl;
+    std::cout << " ThinQRfact MulMat:"<< t4-t3<<" us"<<std::endl;
+  }
+
+  virtual void SolvePrecBlockCG (std::vector<Field> &src, std::vector<Field> &X)
+  {
+    std::cout << GridLogMessage<<"HDCG: mrhs fPrecBlockcg starting"<<std::endl;
+    src[0].Grid()->Barrier();
+    int nrhs = src.size();
+    //    std::vector<RealD> f(nrhs);
+    //    std::vector<RealD> rtzp(nrhs);
+    //    std::vector<RealD> rtz(nrhs);
+    //    std::vector<RealD> a(nrhs);
+    //    std::vector<RealD> d(nrhs);
+    //    std::vector<RealD> b(nrhs);
+    //    std::vector<RealD> rptzp(nrhs);
+
+    ////////////////////////////////////////////
+    //Initial residual computation & set up
+    ////////////////////////////////////////////
+    std::vector<RealD> ssq(nrhs);
+    for(int rhs=0;rhs<nrhs;rhs++){
+      ssq[rhs]=norm2(src[rhs]); assert(ssq[rhs]!=0.0);
+    }      
+
+    ///////////////////////////
+    // Fields -- eliminate duplicates between fPcg and block cg
+    ///////////////////////////
+    std::vector<Field> Mtmp(nrhs,grid);
+    std::vector<Field> tmp(nrhs,grid);
+    std::vector<Field>   Z(nrhs,grid); // Rename Z to R
+    std::vector<Field>  MZ(nrhs,grid); // Rename MZ to Z
+    std::vector<Field>   Q(nrhs,grid); // 
+    std::vector<Field>  MQ(nrhs,grid); // Rename to P
+    std::vector<Field>   D(nrhs,grid);
+    std::vector<Field>  AD(nrhs,grid);
+    
+    /************************************************************************
+     * Preconditioned Block conjugate gradient rQ
+     * Generalise Sebastien Birk Thesis, after Dubrulle 2001.
+     * Introduce preconditioning following Saad Ch9
+     ************************************************************************
+     * Dimensions:
+     *
+     *   X,B etc... ==(Nferm x nrhs)
+     *  Matrix A==(Nferm x Nferm)
+     *  
+     * Nferm = Nspin x Ncolour x Ncomplex x Nlattice_site
+     * QC => Thin QR factorisation (google it)
+     *
+     * R = B-AX
+     * Z = Mi R
+     * QC = Z
+     * D = Q 
+     * for k: 
+     *   R  = AD
+     *   Z  = Mi R
+     *   M  = [D^dag R]^{-1}
+     *   X  = X + D M C
+     *   QS = Q - Z.M
+     *   D  = Q + D S^dag
+     *   C  = S C
+     */
+    Eigen::MatrixXcd m_DZ     = Eigen::MatrixXcd::Identity(nrhs,nrhs);
+    Eigen::MatrixXcd m_M      = Eigen::MatrixXcd::Identity(nrhs,nrhs);
+    Eigen::MatrixXcd m_zz     = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    Eigen::MatrixXcd m_rr     = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    
+    Eigen::MatrixXcd m_C      = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    Eigen::MatrixXcd m_Cinv   = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    Eigen::MatrixXcd m_S      = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    Eigen::MatrixXcd m_Sinv   = Eigen::MatrixXcd::Zero(nrhs,nrhs);
+    
+    Eigen::MatrixXcd m_tmp    = Eigen::MatrixXcd::Identity(nrhs,nrhs);
+    Eigen::MatrixXcd m_tmp1   = Eigen::MatrixXcd::Identity(nrhs,nrhs);
+
+    GridStopWatch HDCGTimer;
+
+    //////////////////////////
+    // x0 = Vstart -- possibly modify guess
+    //////////////////////////
+    Vstart(X,src);
+
+    //////////////////////////
+    // R = B-AX
+    //////////////////////////
+    for(int rhs=0;rhs<nrhs;rhs++){
+      // r0 = b -A x0
+      _FineLinop.HermOp(X[rhs],tmp[rhs]);
+      axpy (Z[rhs], -1.0,tmp[rhs], src[rhs]);    // Computes R=Z=src - A X0
+    }
+
+    //////////////////////////////////
+    // Compute MZ = M1 Z = M1 B - M1 A x0
+    //////////////////////////////////
+    PcgM1(Z,MZ);  
+
+    //////////////////////////////////
+    // QC = Z
+    //////////////////////////////////
+    ThinQRfact (m_zz, m_C, m_Cinv, Q, MQ, Z, MZ);
+
+    //////////////////////////////////
+    // D=MQ
+    //////////////////////////////////
+    for(int b=0;b<nrhs;b++) D[b]=MQ[b]; // LLT rotation of the MZ basis of search dirs
+
+    std::cout << GridLogMessage<<"PrecBlockCGrQ vec computed initial residual and QR fact " <<std::endl;
+
+    ProjectTimer.Reset();
+    PromoteTimer.Reset();
+    DeflateTimer.Reset();
+    CoarseTimer.Reset();
+    SmoothTimer.Reset();
+    FineTimer.Reset();
+    InsertTimer.Reset();
+
+    GridStopWatch M1Timer;
+    GridStopWatch M2Timer;
+    GridStopWatch M3Timer;
+    GridStopWatch LinalgTimer;
+    GridStopWatch InnerProdTimer;
+
+    HDCGTimer.Start();
+
+    std::vector<RealD> rn(nrhs);
+    for (int k=0;k<=MaxIterations;k++){
+
+      ////////////////////
+      // Z  = AD
+      ////////////////////
+      M3Timer.Start();
+      for(int b=0;b<nrhs;b++) _FineLinop.HermOp(D[b], Z[b]);      
+      M3Timer.Stop();
+
+      ////////////////////
+      // MZ  = M1 Z <==== the Multigrid preconditioner
+      ////////////////////
+      M1Timer.Start();
+      PcgM1(Z,MZ);
+      M1Timer.Stop();
+
+      FineTimer.Start();
+      ////////////////////
+      // M  = [D^dag Z]^{-1} = (<Ddag MZ>_M)^{-1} inner prod, generalising Saad derivation of Precon CG
+      ////////////////////
+      InnerProdTimer.Start();
+      _BlockCGLinalg.InnerProductMatrix(m_DZ,D,Z);
+      InnerProdTimer.Stop();
+      m_M       = m_DZ.inverse();
+
+      ///////////////////////////
+      // X  = X + D MC
+      ///////////////////////////
+      m_tmp     = m_M * m_C;
+      LinalgTimer.Start();
+      _BlockCGLinalg.MaddMatrix(X,m_tmp, D,X);     // D are the search directions and X takes the updates 
+      LinalgTimer.Stop();
+
+      ///////////////////////////
+      // QS = Q - M Z
+      // (MQ) S = MQ - M (M1Z)
+      ///////////////////////////
+      LinalgTimer.Start();
+      _BlockCGLinalg.MaddMatrix(tmp ,m_M, Z, Q,-1.0);
+      _BlockCGLinalg.MaddMatrix(Mtmp,m_M,MZ,MQ,-1.0);
+      ThinQRfact (m_zz, m_S, m_Sinv, Q, MQ, tmp, Mtmp);
+      LinalgTimer.Stop();
+
+      ////////////////////////////
+      // D  = MQ + D S^dag
+      ////////////////////////////
+      m_tmp = m_S.adjoint();
+      LinalgTimer.Start();
+      _BlockCGLinalg.MaddMatrix(D,m_tmp,D,MQ);
+      LinalgTimer.Stop();
+
+      ////////////////////////////
+      // C  = S C
+      ////////////////////////////
+      m_C = m_S*m_C;
+      
+      ////////////////////////////
+      // convergence monitor
+      ////////////////////////////
+      m_rr = m_C.adjoint() * m_C;
+      
+      FineTimer.Stop();
+
+      RealD max_resid=0;
+      RealD rrsum=0;
+      RealD sssum=0;
+      RealD rr;
+
+      for(int b=0;b<nrhs;b++) {
+	rrsum+=real(m_rr(b,b));
+	sssum+=ssq[b];
+	rr = real(m_rr(b,b))/ssq[b];
+	if ( rr > max_resid ) max_resid = rr;
+      }
+      std::cout << GridLogMessage <<
+	  "\t Prec BlockCGrQ Iteration "<<k<<" ave resid "<< std::sqrt(rrsum/sssum) << " max "<< std::sqrt(max_resid) <<std::endl;
+
+
+      if ( max_resid < Tolerance*Tolerance ) { 
+
+	HDCGTimer.Stop();
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ converged in "<<k<<" iterations and "<<HDCGTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Linalg  "<<LinalgTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : fine H  "<<M3Timer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : prec M1 "<<M1Timer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"**** M1 breakdown:"<<std::endl;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Project "<<ProjectTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Promote "<<PromoteTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Deflate "<<DeflateTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Coarse  "<<CoarseTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Fine    "<<FineTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Smooth  "<<SmoothTimer.Elapsed()<<std::endl;;
+	std::cout<<GridLogMessage<<"HDCG: mrhs PrecBlockCGrQ : Insert  "<<InsertTimer.Elapsed()<<std::endl;;
+
+	for(int rhs=0;rhs<nrhs;rhs++){
+
+	  _FineLinop.HermOp(X[rhs],tmp[rhs]);			  
+
+	  Field mytmp(grid);
+	  axpy(mytmp,-1.0,src[rhs],tmp[rhs]);
+      
+	  RealD  xnorm   = sqrt(norm2(X[rhs]));
+	  RealD  srcnorm = sqrt(norm2(src[rhs]));
+	  RealD  tmpnorm = sqrt(norm2(mytmp));
+	  RealD  true_residual = tmpnorm/srcnorm;
+	  std::cout<<GridLogMessage
+		   <<"HDCG: true residual ["<<rhs<<"] is "<<true_residual
+		   <<" solution "<<xnorm
+		   <<" source "<<srcnorm
+		   <<std::endl;
+	}
+	return;
+      }
+      
+    }
+    HDCGTimer.Stop();
+    std::cout<<GridLogMessage<<"HDCG: PrecBlockCGrQ not converged "<<HDCGTimer.Elapsed()<<std::endl;
+    assert(0);
+  }
+
+  virtual void SolveSingleSystem (std::vector<Field> &src, std::vector<Field> &x)
   {
     std::cout << GridLogMessage<<"HDCG: mrhs fPcg starting"<<std::endl;
     src[0].Grid()->Barrier();
@@ -361,15 +657,23 @@ public:
     CoarseField PleftProjMrhs(this->coarsegridmrhs);
     CoarseField PleftMss_projMrhs(this->coarsegridmrhs);
 
+#undef SMOOTHER_BLOCK_SOLVE
+#if SMOOTHER_BLOCK_SOLVE
+    this->SmoothTimer.Start();
+    this->_Smoother(in,Min);
+    this->SmoothTimer.Stop();
+#else
     for(int rhs=0;rhs<nrhs;rhs++) {
-
       this->SmoothTimer.Start();
       this->_Smoother(in[rhs],Min[rhs]);
       this->SmoothTimer.Stop();
-
+    }
+#endif
+    
+    for(int rhs=0;rhs<nrhs;rhs++) {
+      
       this->FineTimer.Start();
       this->_FineLinop.HermOp(Min[rhs],out[rhs]);
-
       axpy(tmp[rhs],-1.0,out[rhs],in[rhs]);          // resid  = in - A Min
       this->FineTimer.Stop();
 
@@ -407,7 +711,7 @@ public:
     this->FineTimer.Stop();
   }
 };
-  
+
 
 NAMESPACE_END(Grid);
 

Grid/algorithms/iterative/BlockConjugateGradient.h

@@ -31,6 +31,58 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
+template<class Field>
+void InnerProductMatrix(Eigen::MatrixXcd &m , const std::vector<Field> &X, const std::vector<Field> &Y){
+  typedef typename Field::scalar_type scomplex;
+  int Nblock = X.size();
+  for(int b=0;b<Nblock;b++){
+  for(int bp=0;bp<Nblock;bp++) {
+    m(b,bp) = innerProduct(X[b],Y[bp]);  
+  }}
+}
+template<class Field>
+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
+  //
+  //Could pack "X" and "AP" into a Nblock x Volume dense array.
+  // AP(Nrhs x vol) = Y(Nrhs x vol) + scale * m(nrhs x nrhs) * X(nrhs*vol)
+  typedef typename Field::scalar_type scomplex;
+  int Nblock = AP.size();
+  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] +scomplex(scale*m(bp,b))*X[bp]; 
+    }
+  }
+  for(int b=0;b<Nblock;b++){
+    AP[b] = tmp[b];
+  }
+}
+template<class Field>
+void MulMatrix(std::vector<Field> &AP, Eigen::MatrixXcd &m , const std::vector<Field> &X){
+  // Should make this cache friendly with site outermost, parallel_for
+  typedef typename Field::scalar_type scomplex;
+  int Nblock = AP.size();
+  for(int b=0;b<Nblock;b++){
+    AP[b] = Zero();
+    for(int bp=0;bp<Nblock;bp++) {
+      AP[b] += scomplex(m(bp,b))*X[bp]; 
+    }
+  }
+}
+template<class Field>
+double normv(const std::vector<Field> &P){
+  int Nblock = P.size();
+  double nn = 0.0;
+  for(int b=0;b<Nblock;b++) {
+    nn+=norm2(P[b]);
+  }
+  return nn;
+}
+
+
 enum BlockCGtype { BlockCG, BlockCGrQ, CGmultiRHS, BlockCGVec, BlockCGrQVec };
 
 //////////////////////////////////////////////////////////////////////////
@@ -87,10 +139,19 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr,
   sliceInnerProductMatrix(m_rr,R,R,Orthog);
 
   // Force manifest hermitian to avoid rounding related
+  /*
+  int rank=m_rr.rows();
+  for(int r=0;r<rank;r++){
+  for(int s=0;s<rank;s++){
+    std::cout << "QR m_rr["<<r<<","<<s<<"] "<<m_rr(r,s)<<std::endl;
+  }}
+  */
   m_rr = 0.5*(m_rr+m_rr.adjoint());
 
   Eigen::MatrixXcd L    = m_rr.llt().matrixL(); 
 
+//  ComplexD det = L.determinant();
+//  std::cout << " Det m_rr "<<det<<std::endl;
   C    = L.adjoint();
   Cinv = C.inverse();
   ////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -110,11 +171,20 @@ void ThinQRfact (Eigen::MatrixXcd &m_rr,
 		 const std::vector<Field> & R)
 {
   InnerProductMatrix(m_rr,R,R);
-
+  /*
+  int rank=m_rr.rows();
+  for(int r=0;r<rank;r++){
+  for(int s=0;s<rank;s++){
+    std::cout << "QRvec m_rr["<<r<<","<<s<<"] "<<m_rr(r,s)<<std::endl;
+  }}
+  */
   m_rr = 0.5*(m_rr+m_rr.adjoint());
 
   Eigen::MatrixXcd L    = m_rr.llt().matrixL(); 
 
+  //  ComplexD det = L.determinant();
+  //  std::cout << " Det m_rr "<<det<<std::endl;
+
   C    = L.adjoint();
   Cinv = C.inverse();
 
@@ -186,6 +256,7 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   sliceNorm(ssq,B,Orthog);
   RealD sssum=0;
   for(int b=0;b<Nblock;b++) sssum+=ssq[b];
+  for(int b=0;b<Nblock;b++) std::cout << "src["<<b<<"]" << ssq[b] <<std::endl;
 
   sliceNorm(residuals,B,Orthog);
   for(int b=0;b<Nblock;b++){ assert(std::isnan(residuals[b])==0); }
@@ -221,6 +292,9 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   Linop.HermOp(X, AD);
   tmp = B - AD;  
 
+  sliceNorm(residuals,tmp,Orthog);
+  for(int b=0;b<Nblock;b++) std::cout << "res["<<b<<"]" << residuals[b] <<std::endl;
+  
   ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);
   D=Q;
 
@@ -236,6 +310,8 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
   GridStopWatch SolverTimer;
   SolverTimer.Start();
 
+  RealD max_resid=0;
+
   int k;
   for (k = 1; k <= MaxIterations; k++) {
 
@@ -280,7 +356,7 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
      */
     m_rr = m_C.adjoint() * m_C;
 
-    RealD max_resid=0;
+    max_resid=0;
     RealD rrsum=0;
     RealD rr;
 
@@ -322,7 +398,9 @@ void BlockCGrQsolve(LinearOperatorBase<Field> &Linop, const Field &B, Field &X)
     }
 
   }
-  std::cout << GridLogMessage << "BlockConjugateGradient(rQ) did NOT converge" << std::endl;
+
+  std::cout << GridLogMessage << "BlockConjugateGradient(rQ) did NOT converge "<<k<<" / "<<MaxIterations
+	    <<" residual "<< std::sqrt(max_resid)<< std::endl;
 
   if (ErrorOnNoConverge) assert(0);
   IterationsToComplete = k;
@@ -466,43 +544,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] + scomplex(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] += scomplex(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:
 //--------------------------
@@ -549,6 +590,7 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
 
   RealD sssum=0;
   for(int b=0;b<Nblock;b++){ ssq[b] = norm2(B[b]);}
+  for(int b=0;b<Nblock;b++){ std::cout << "ssq["<<b<<"] "<<ssq[b]<<std::endl;}
   for(int b=0;b<Nblock;b++) sssum+=ssq[b];
 
   for(int b=0;b<Nblock;b++){ residuals[b] = norm2(B[b]);}
@@ -585,6 +627,7 @@ void BlockCGrQsolveVec(LinearOperatorBase<Field> &Linop, const std::vector<Field
   for(int b=0;b<Nblock;b++) {
     Linop.HermOp(X[b], AD[b]);
     tmp[b] = B[b] - AD[b];  
+    std::cout << "r0["<<b<<"] "<<norm2(tmp[b])<<std::endl;
   }
 
   ThinQRfact (m_rr, m_C, m_Cinv, Q, tmp);

Grid/algorithms/iterative/ConjugateGradient.h

@@ -38,12 +38,13 @@ NAMESPACE_BEGIN(Grid);
 // single input vec, single output vec.
 /////////////////////////////////////////////////////////////
 
+
 template <class Field>
 class ConjugateGradient : public OperatorFunction<Field> {
 public:
 
   using OperatorFunction<Field>::operator();
-
+  
   bool ErrorOnNoConverge;  // throw an assert when the CG fails to converge.
                            // Defaults true.
   RealD Tolerance;
@@ -57,10 +58,22 @@ public:
       ErrorOnNoConverge(err_on_no_conv)
   {};
 
-  void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
+  virtual void LogIteration(int k,RealD a,RealD b){
+    //    std::cout << "ConjugageGradient::LogIteration() "<<std::endl;
+  };
+  virtual void LogBegin(void){
+    std::cout << "ConjugageGradient::LogBegin() "<<std::endl;
+  };
 
-    GRID_TRACE("ConjugateGradient");
+    void operator()(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi) {
+
+      this->LogBegin();
+
+      GRID_TRACE("ConjugateGradient");
     GridStopWatch PreambleTimer;
+    GridStopWatch ConstructTimer;
+    GridStopWatch NormTimer;
+    GridStopWatch AssignTimer;
     PreambleTimer.Start();
     psi.Checkerboard() = src.Checkerboard();
 
@@ -70,14 +83,19 @@ public:
     //RealD b_pred;
 
     // Was doing copies
-    Field p(src.Grid());
+    ConstructTimer.Start();
+    Field p  (src.Grid());
     Field mmp(src.Grid());
-    Field r(src.Grid());
+    Field r  (src.Grid());
+    ConstructTimer.Stop();
 
     // Initial residual computation & set up
+    NormTimer.Start();
     ssq = norm2(src);
     RealD guess = norm2(psi);
+    NormTimer.Stop();
     assert(std::isnan(guess) == 0);
+    AssignTimer.Start();
     if ( guess == 0.0 ) {
       r = src;
       p = r;
@@ -89,6 +107,7 @@ public:
       a = norm2(p);
     }
     cp = a;
+    AssignTimer.Stop();
 
     // Handle trivial case of zero src
     if (ssq == 0.){
@@ -164,6 +183,7 @@ public:
       }
       LinearCombTimer.Stop();
       LinalgTimer.Stop();
+      LogIteration(k,a,b);
 
       IterationTimer.Stop();
       if ( (k % 500) == 0 ) {
@@ -220,6 +240,9 @@ public:
     	      <<" residual "<< std::sqrt(cp / ssq)<< std::endl;
     SolverTimer.Stop();
     std::cout << GridLogMessage << "\tPreamble   " << PreambleTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tConstruct  " << ConstructTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tNorm       " << NormTimer.Elapsed() <<std::endl;
+    std::cout << GridLogMessage << "\tAssign     " << AssignTimer.Elapsed() <<std::endl;
     std::cout << GridLogMessage << "\tSolver     " << SolverTimer.Elapsed() <<std::endl;
     std::cout << GridLogMessage << "Solver breakdown "<<std::endl;
     std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
@@ -233,5 +256,118 @@ public:
 
   }
 };
+
+
+template <class Field>
+class ConjugateGradientPolynomial : public ConjugateGradient<Field> {
+public:
+  // Optionally record the CG polynomial
+  std::vector<double> ak;
+  std::vector<double> bk;
+  std::vector<double> poly_p;
+  std::vector<double> poly_r;
+  std::vector<double> poly_Ap;
+  std::vector<double> polynomial;
+
+public:
+  ConjugateGradientPolynomial(RealD tol, Integer maxit, bool err_on_no_conv = true)
+    : ConjugateGradient<Field>(tol,maxit,err_on_no_conv)
+  { };
+  void PolyHermOp(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi)
+  {
+    Field tmp(src.Grid());
+    Field AtoN(src.Grid());
+    AtoN = src;
+    psi=AtoN*polynomial[0];
+    for(int n=1;n<polynomial.size();n++){
+      tmp = AtoN;
+      Linop.HermOp(tmp,AtoN);
+      psi = psi + polynomial[n]*AtoN;
+    }
+  }
+  void CGsequenceHermOp(LinearOperatorBase<Field> &Linop, const Field &src, Field &x)
+  {
+    Field Ap(src.Grid());
+    Field r(src.Grid());
+    Field p(src.Grid());
+    p=src;
+    r=src;
+    x=Zero();
+    x.Checkerboard()=src.Checkerboard();
+    for(int k=0;k<ak.size();k++){
+      x = x + ak[k]*p;
+      Linop.HermOp(p,Ap);
+      r = r - ak[k] * Ap;
+      p = r + bk[k] * p;
+    }
+  }
+  void Solve(LinearOperatorBase<Field> &Linop, const Field &src, Field &psi)
+  {
+    psi=Zero();
+    this->operator ()(Linop,src,psi);
+  }
+  virtual void LogBegin(void)
+  {
+    std::cout << "ConjugageGradientPolynomial::LogBegin() "<<std::endl;
+    ak.resize(0);
+    bk.resize(0);
+    polynomial.resize(0);
+    poly_Ap.resize(0);
+    poly_Ap.resize(0);
+    poly_p.resize(1);
+    poly_r.resize(1);
+    poly_p[0]=1.0;
+    poly_r[0]=1.0;
+  };
+  virtual void LogIteration(int k,RealD a,RealD b)
+  {
+    // With zero guess,
+    // p = r = src
+    //
+    // iterate:
+    //   x =  x + a p
+    //   r =  r - a A p
+    //   p =  r + b p
+    //
+    // [0]
+    // r = x
+    // p = x
+    // Ap=0
+    //
+    // [1]
+    // Ap = A x + 0  ==> shift poly P right by 1 and add 0.
+    // x  = x + a p  ==> add polynomials term by term 
+    // r  = r - a A p  ==> add polynomials term by term
+    // p  = r + b p  ==> add polynomials term by term
+    //
+    std::cout << "ConjugageGradientPolynomial::LogIteration() "<<k<<std::endl;
+    ak.push_back(a);
+    bk.push_back(b);
+    //  Ap= right_shift(p)
+    poly_Ap.resize(k+1);
+    poly_Ap[0]=0.0;
+    for(int i=0;i<k;i++){
+      poly_Ap[i+1]=poly_p[i];
+    }
+
+    //  x = x + a p
+    polynomial.resize(k);
+    polynomial[k-1]=0.0;
+    for(int i=0;i<k;i++){
+      polynomial[i] = polynomial[i] + a * poly_p[i];
+    }
+    
+    //  r = r - a Ap
+    //  p = r + b p
+    poly_r.resize(k+1);
+    poly_p.resize(k+1);
+    poly_r[k] = poly_p[k] = 0.0;
+    for(int i=0;i<k+1;i++){
+      poly_r[i] = poly_r[i] - a * poly_Ap[i];
+      poly_p[i] = poly_r[i] + b * poly_p[i];
+    }
+  }
+};
+
 NAMESPACE_END(Grid);
 #endif

Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h

@@ -143,7 +143,7 @@ public:
       ip = innerProduct(evec[j],w); 
       if(if_print) 
       if( norm(ip)/norm2(w) > 1e-14)
-      Glog<<"orthogonalize before: "<<j<<" of "<<k<<" "<< ip <<std::endl;
+	Glog<<"orthogonalize before: "<<j<<" of "<<k<<" "<< ip <<std::endl;
       w = w - ip * evec[j];
       if(if_print) {
         ip = innerProduct(evec[j],w); 
@@ -279,16 +279,16 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
       diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
       _sort.push(eval2,Nm);
-      Glog << "#Ritz value before shift: "<< std::endl;
+      //      Glog << "#Ritz value before shift: "<< std::endl;
       for(int i=0; i<Nm; ++i){
-	std::cout.precision(13);
-	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-	std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+	//	std::cout.precision(13);
+	//	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	//	std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
       }
       
       //----------------------------------------------------------------------
       if ( Nm>Nk ) {
-        Glog <<" #Apply shifted QR transformations "<<std::endl;
+	//        Glog <<" #Apply shifted QR transformations "<<std::endl;
         //int k2 = Nk+Nu;
         int k2 = Nk;
       
@@ -326,11 +326,11 @@ public:
         Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
         diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
         _sort.push(eval2,Nk);
-	Glog << "#Ritz value after shift: "<< std::endl;
+	//	Glog << "#Ritz value after shift: "<< std::endl;
         for(int i=0; i<Nk; ++i){
-	  //          std::cout.precision(13);
-	  //          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
-	  //          std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
+	  //	  std::cout.precision(13);
+	  //	  std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	  //	  std::cout << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
         }
       }
       //----------------------------------------------------------------------
@@ -644,7 +644,7 @@ private:
       //      for (int u=0; u<mrhs; ++u) Glog << " out["<<u<<"] = "<<norm2(out[u])<<std::endl;
       k_start +=mrhs;
     }
-    Glog << "LinAlg "<< std::endl;
+    //    Glog << "LinAlg "<< std::endl;
     
     if (b>0) {
       for (int u=0; u<Nu; ++u) {
@@ -678,7 +678,7 @@ private:
       }
       w_copy[u] = w[u];
     }
-    Glog << "LinAlg done"<< std::endl;
+    //    Glog << "LinAlg done"<< std::endl;
     
     // In block version, the steps 6 and 7 in Lanczos construction is
     // replaced by the QR decomposition of new basis block.
@@ -691,15 +691,15 @@ private:
     }
 
     // re-orthogonalization for numerical stability
-    Glog << "Gram Schmidt"<< std::endl;
+    //    Glog << "Gram Schmidt"<< std::endl;
     orthogonalize(w,Nu,evec,R);
     // QR part
     for (int u=1; u<Nu; ++u) {
       orthogonalize(w[u],w,u);
     }
-    Glog << "Gram Schmidt done "<< std::endl;
+    //    Glog << "Gram Schmidt done "<< std::endl;
     
-    Glog << "LinAlg "<< std::endl;
+    //    Glog << "LinAlg "<< std::endl;
     for (int u=0; u<Nu; ++u) {
       //for (int v=0; v<Nu; ++v) {
       for (int v=u; v<Nu; ++v) {
@@ -716,7 +716,7 @@ private:
 	//        Glog <<" In block "<< b << "," <<" beta[" << u << "," << k-L << "] = " << lme[u][k] << std::endl;
       }
     }
-    Glog << "LinAlg done "<< std::endl;
+    //    Glog << "LinAlg done "<< std::endl;
 
     if (b < Nm/Nu-1) {
       for (int u=0; u<Nu; ++u) {
@@ -935,7 +935,7 @@ if (1){
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
+    //    Glog << "unpackHermitBlockTriDiagMatToEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     M = Eigen::MatrixXcd::Zero(Nk,Nk);
@@ -953,7 +953,7 @@ if (1){
         M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
       }
     }
-    Glog << "unpackHermitBlockTriDiagMatToEigen() end" << std::endl; 
+    //    Glog << "unpackHermitBlockTriDiagMatToEigen() end" << std::endl; 
   }
  
 
@@ -963,7 +963,7 @@ if (1){
          int Nu, int Nb, int Nk, int Nm,
          Eigen::MatrixXcd& M)
   {
-    Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
+    //    Glog << "packHermitBlockTriDiagMatfromEigen() begin" << '\n'; 
     assert( Nk%Nu == 0 && Nm%Nu == 0 );
     assert( Nk <= Nm );
     
@@ -979,7 +979,7 @@ if (1){
         lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
       }
     }
-    Glog << "packHermitBlockTriDiagMatfromEigen() end" <<std::endl; 
+    //    Glog << "packHermitBlockTriDiagMatfromEigen() end" <<std::endl; 
   }
 
 
@@ -988,7 +988,7 @@ if (1){
 		            RealD Dsh,
 		            Eigen::MatrixXcd& Qprod)
   {
-    Glog << "shiftedQRDecompEigen() begin" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() begin" << '\n'; 
     Eigen::MatrixXcd Q = Eigen::MatrixXcd::Zero(Nm,Nm);
     Eigen::MatrixXcd R = Eigen::MatrixXcd::Zero(Nm,Nm);
     Eigen::MatrixXcd Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
@@ -1004,7 +1004,7 @@ if (1){
                         // lower triangular part used to represent series
                         // of Q sequence.
 
-    Glog << "shiftedQRDecompEigen() Housholder & QR" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Housholder & QR" << '\n'; 
     // equivalent operation of Qprod *= Q
     //M = Eigen::MatrixXcd::Zero(Nm,Nm);
     
@@ -1025,7 +1025,7 @@ if (1){
     
     Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
 
-    Glog << "shiftedQRDecompEigen() Mtmp create" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Mtmp create" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=0; j<Nm-(Nu+1); ++j) {
         for (int k=0; k<Nu+1+j; ++k) {
@@ -1033,7 +1033,7 @@ if (1){
         }
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop1" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Mtmp loop1" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=Nm-(Nu+1); j<Nm; ++j) {
         for (int k=0; k<Nm; ++k) {
@@ -1041,7 +1041,7 @@ if (1){
         }
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop2" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Mtmp loop2" << '\n'; 
     
     //static int ntimes = 2;
     //for (int j=0; j<Nm-(ntimes*Nu); ++j) {
@@ -1067,13 +1067,13 @@ if (1){
         Mtmp(j,i) = conj(Mtmp(i,j));
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop3" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Mtmp loop3" << '\n'; 
 
     for (int i=0; i<Nm; ++i) {
       Mtmp(i,i) = real(Mtmp(i,i)) + Dsh;
     }
     
-    Glog << "shiftedQRDecompEigen() Mtmp loop4" << '\n'; 
+    //    Glog << "shiftedQRDecompEigen() Mtmp loop4" << '\n'; 
     M = Mtmp;
 
     //M = Q.adjoint()*(M*Q);
@@ -1085,7 +1085,7 @@ if (1){
     //  }
     //}
     
-    Glog << "shiftedQRDecompEigen() end" <<std::endl; 
+    //    Glog << "shiftedQRDecompEigen() end" <<std::endl; 
   }
 
   void exampleQRDecompEigen(void)

Grid/algorithms/iterative/NormalEquations.h

@@ -60,6 +60,32 @@ public:
   }     
 };
 
+template<class Field> class NormalResidual : public LinearFunction<Field>{
+private:
+  SparseMatrixBase<Field> & _Matrix;
+  OperatorFunction<Field> & _HermitianSolver;
+  LinearFunction<Field>   & _Guess;
+public:
+
+  /////////////////////////////////////////////////////
+  // Wrap the usual normal equations trick
+  /////////////////////////////////////////////////////
+ NormalResidual(SparseMatrixBase<Field> &Matrix, OperatorFunction<Field> &HermitianSolver,
+		 LinearFunction<Field> &Guess) 
+   :  _Matrix(Matrix), _HermitianSolver(HermitianSolver), _Guess(Guess) {}; 
+
+  void operator() (const Field &in, Field &out){
+ 
+    Field res(in.Grid());
+    Field tmp(in.Grid());
+
+    MMdagLinearOperator<SparseMatrixBase<Field>,Field> MMdagOp(_Matrix);
+    _Guess(in,res);
+    _HermitianSolver(MMdagOp,in,res);  // M Mdag res = in ;
+    _Matrix.Mdag(res,out);             // out = Mdag res
+  }     
+};
+
 template<class Field> class HPDSolver : public LinearFunction<Field> {
 private:
   LinearOperatorBase<Field> & _Matrix;

Grid/algorithms/iterative/PowerMethod.h

@@ -20,7 +20,7 @@ template<class Field> class PowerMethod
     RealD evalMaxApprox = 0.0; 
     auto src_n = src; 
     auto tmp = src; 
-    const int _MAX_ITER_EST_ = 100; 
+    const int _MAX_ITER_EST_ = 200; 
 
     for (int i=0;i<_MAX_ITER_EST_;i++) { 
       
@@ -30,18 +30,17 @@ template<class Field> class PowerMethod
       RealD vden = norm2(src_n); 
       RealD na = vnum/vden; 
 
-      std::cout << GridLogIterative << "PowerMethod: Current approximation of largest eigenvalue " << na << std::endl;
+      std::cout << GridLogMessage << "PowerMethod: Current approximation of largest eigenvalue " << na << std::endl;
       
-      if ( (fabs(evalMaxApprox/na - 1.0) < 0.001) || (i==_MAX_ITER_EST_-1) ) { 
- 	evalMaxApprox = na; 
-	std::cout << GridLogMessage << " Approximation of largest eigenvalue: " << evalMaxApprox << std::endl;
- 	return evalMaxApprox; 
-      } 
+      //      if ( (fabs(evalMaxApprox/na - 1.0) < 0.0001) || (i==_MAX_ITER_EST_-1) ) { 
+	// 	evalMaxApprox = na; 
+	// 	return evalMaxApprox; 
+      //      } 
       evalMaxApprox = na; 
       src_n = tmp;
     }
-    assert(0);
-    return 0;
+    std::cout << GridLogMessage << " Approximation of largest eigenvalue: " << evalMaxApprox << std::endl;
+    return evalMaxApprox;
   }
 };
 }

Grid/algorithms/iterative/PowerSpectrum.h

@@ -0,0 +1,76 @@
+#pragma once
+namespace Grid {
+
+class Band
+{
+  RealD lo, hi;
+public:
+  Band(RealD _lo,RealD _hi)
+  {
+    lo=_lo;
+    hi=_hi;
+  }
+  RealD operator() (RealD x){
+    if ( x>lo && x<hi ){
+      return 1.0;
+    } else {
+      return 0.0;
+    }
+  }
+};
+
+class PowerSpectrum
+{ 
+ public: 
+
+  template<typename T>  static RealD normalise(T& v) 
+  {
+    RealD nn = norm2(v);
+    nn = sqrt(nn);
+    v = v * (1.0/nn);
+    return nn;
+  }
+
+  std::vector<RealD> ranges;
+  std::vector<int> order;
+  
+  PowerSpectrum(  std::vector<RealD> &bins, std::vector<int> &_order ) : ranges(bins), order(_order)  { };
+
+  template<class Field>
+  RealD operator()(LinearOperatorBase<Field> &HermOp, const Field &src) 
+  { 
+    GridBase *grid = src.Grid(); 
+    int N=ranges.size();
+    RealD hi = ranges[N-1];
+
+    RealD lo_band = 0.0;
+    RealD hi_band;
+    RealD nn=norm2(src);
+    RealD ss=0.0;
+
+    Field tmp = src;
+
+    for(int b=0;b<N;b++){
+      hi_band = ranges[b];
+      Band Notch(lo_band,hi_band);
+      
+      Chebyshev<Field> polynomial;
+      polynomial.Init(0.0,hi,order[b],Notch);
+      polynomial.JacksonSmooth();
+
+      polynomial(HermOp,src,tmp) ;
+
+      RealD p=norm2(tmp);
+      ss=ss+p;
+      std::cout << GridLogMessage << " PowerSpectrum Band["<<lo_band<<","<<hi_band<<"] power "<<norm2(tmp)/nn<<std::endl;
+      
+      lo_band=hi_band;
+    }
+    std::cout << GridLogMessage << " PowerSpectrum total power "<<ss/nn<<std::endl;
+    std::cout << GridLogMessage << " PowerSpectrum total power (unnormalised) "<<nn<<std::endl;
+
+    return 0;
+  };
+};
+  
+}

Grid/cartesian/Cartesian_base.h

@@ -82,6 +82,7 @@ public:
   bool _isCheckerBoarded; 
   int        LocallyPeriodic;
   Coordinate _checker_dim_mask;
+  int              _checker_dim;
 
 public:
 
@@ -89,7 +90,7 @@ public:
   // Checkerboarding interface is virtual and overridden by 
   // GridCartesian / GridRedBlackCartesian
   ////////////////////////////////////////////////////////////////
-  virtual int CheckerBoarded(int dim)=0;
+  virtual int CheckerBoarded(int dim) =0;
   virtual int CheckerBoard(const Coordinate &site)=0;
   virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;
   virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0;

Grid/cartesian/Cartesian_full.h

@@ -38,7 +38,7 @@ class GridCartesian: public GridBase {
 
 public:
   int dummy;
-  Coordinate _checker_dim_mask;
+  //  Coordinate _checker_dim_mask;
   virtual int  CheckerBoardFromOindexTable (int Oindex) {
     return 0;
   }
@@ -46,7 +46,7 @@ public:
   {
     return 0;
   }
-  virtual int CheckerBoarded(int dim){
+  virtual int CheckerBoarded(int dim) {
     return 0;
   }
   virtual int CheckerBoard(const Coordinate &site){
@@ -106,6 +106,7 @@ public:
     _rdimensions.resize(_ndimension);
     _simd_layout.resize(_ndimension);
     _checker_dim_mask.resize(_ndimension);;
+    _checker_dim = -1;
     _lstart.resize(_ndimension);
     _lend.resize(_ndimension);
 

Grid/cartesian/Cartesian_red_black.h

@@ -57,9 +57,10 @@ class GridRedBlackCartesian : public GridBase
 {
 public:
   //  Coordinate _checker_dim_mask;
-  int              _checker_dim;
+  //  int              _checker_dim;
   std::vector<int> _checker_board;
 
+  virtual int isCheckerBoarded(void) const { return 1; };
   virtual int CheckerBoarded(int dim){
     if( dim==_checker_dim) return 1;
     else return 0;
@@ -147,7 +148,7 @@ public:
   {
     Init(base->_fdimensions,base->_simd_layout,base->_processors,checker_dim_mask,checker_dim)  ;
   }
-
+  
   virtual ~GridRedBlackCartesian() = default;
 
   void Init(const Coordinate &dimensions,

Grid/lattice/Lattice_base.h

@@ -236,17 +236,17 @@ public:
   template<class sobj> inline Lattice<vobj> & operator = (const sobj & r){
     vobj vtmp;
     vtmp = r;
-#if 1
-    auto me  = View(CpuWrite);
-    thread_for(ss,me.size(),{
-       me[ss]= r;
-      });
-#else    
+#if defined(GRID_HIP) || defined(GRID_CUDA) || defined (GRID_SYCL)
     auto me  = View(AcceleratorWrite);
     accelerator_for(ss,me.size(),vobj::Nsimd(),{
 	auto stmp=coalescedRead(vtmp);
 	coalescedWrite(me[ss],stmp);
     });
+#else    
+    auto me  = View(CpuWrite);
+    thread_for(ss,me.size(),{
+       me[ss]= r;
+      });
 #endif    
     me.ViewClose();
     return *this;

Grid/lattice/Lattice_peekpoke.h

@@ -165,7 +165,7 @@ inline void peekLocalSite(sobj &s,const LatticeView<vobj> &l,Coordinate &site)
 
   int Nsimd = grid->Nsimd();
 
-  assert( l.Checkerboard()== grid->CheckerBoard(site));
+  //  assert( l.Checkerboard()== grid->CheckerBoard(site));
   assert( sizeof(sobj)*Nsimd == sizeof(vobj));
 
   static const int words=sizeof(vobj)/sizeof(vector_type);
@@ -179,7 +179,7 @@ inline void peekLocalSite(sobj &s,const LatticeView<vobj> &l,Coordinate &site)
   for(int w=0;w<words;w++){
     pt[w] = getlane(vp[w],idx);
   }
-      
+  //  std::cout << "peekLocalSite "<<site<<" "<<odx<<","<<idx<<" "<<s<<std::endl;
   return;
 };
 template<class vobj,class sobj>
@@ -202,7 +202,7 @@ inline void pokeLocalSite(const sobj &s,LatticeView<vobj> &l,Coordinate &site)
 
   int Nsimd = grid->Nsimd();
 
-  assert( l.Checkerboard()== grid->CheckerBoard(site));
+  //  assert( l.Checkerboard()== grid->CheckerBoard(site));
   assert( sizeof(sobj)*Nsimd == sizeof(vobj));
 
   static const int words=sizeof(vobj)/sizeof(vector_type);

Grid/lattice/Lattice_reduction.h

@@ -519,7 +519,20 @@ sliceSum(const Lattice<vobj> &Data,int orthogdim)
   return result;
 }
 
+/*
+Reimplement
 
+1)
+template<class vobj>
+static void sliceMaddMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<vobj> &X,const Lattice<vobj> &Y,int Orthog,RealD scale=1.0) 
+
+2)
+template<class vobj>
+static void sliceInnerProductMatrix(  Eigen::MatrixXcd &mat, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog) 
+
+3)
+-- Make Slice Mul Matrix call sliceMaddMatrix
+ */
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {
@@ -670,203 +683,96 @@ static void sliceMaddVector(Lattice<vobj> &R,std::vector<RealD> &a,const Lattice
   }
 };
 
-/*
 inline GridBase         *makeSubSliceGrid(const GridBase *BlockSolverGrid,int Orthog)
 {
   int NN    = BlockSolverGrid->_ndimension;
   int nsimd = BlockSolverGrid->Nsimd();
   
-  std::vector<int> latt_phys(0);
-  std::vector<int> simd_phys(0);
-  std::vector<int>  mpi_phys(0);
-  
+  std::vector<int> latt_phys(NN-1);
+  Coordinate simd_phys;
+  std::vector<int>  mpi_phys(NN-1);
+  Coordinate checker_dim_mask(NN-1);
+  int checker_dim=-1;
+
+  int dd;
   for(int d=0;d<NN;d++){
     if( d!=Orthog ) { 
-      latt_phys.push_back(BlockSolverGrid->_fdimensions[d]);
-      simd_phys.push_back(BlockSolverGrid->_simd_layout[d]);
-      mpi_phys.push_back(BlockSolverGrid->_processors[d]);
+      latt_phys[dd]=BlockSolverGrid->_fdimensions[d];
+      mpi_phys[dd] =BlockSolverGrid->_processors[d];
+      checker_dim_mask[dd] = BlockSolverGrid->_checker_dim_mask[d];
+      if ( d == BlockSolverGrid->_checker_dim ) checker_dim = dd;
+      dd++;
     }
   }
-  return (GridBase *)new GridCartesian(latt_phys,simd_phys,mpi_phys); 
+  simd_phys=GridDefaultSimd(latt_phys.size(),nsimd);
+  GridCartesian *tmp         = new GridCartesian(latt_phys,simd_phys,mpi_phys);
+  if(BlockSolverGrid->_isCheckerBoarded) {
+    GridRedBlackCartesian *ret = new GridRedBlackCartesian(tmp,checker_dim_mask,checker_dim);
+    delete tmp;
+    return (GridBase *) ret;
+  } else { 
+    return (GridBase *) tmp;
+  }
 }
-*/
 
 template<class vobj>
 static void sliceMaddMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<vobj> &X,const Lattice<vobj> &Y,int Orthog,RealD scale=1.0) 
 {    
+  GridBase *FullGrid = X.Grid();
+  GridBase *SliceGrid = makeSubSliceGrid(FullGrid,Orthog);
+
+  Lattice<vobj> Ys(SliceGrid);
+  Lattice<vobj> Rs(SliceGrid);
+  Lattice<vobj> Xs(SliceGrid);
+  Lattice<vobj> RR(FullGrid);
+
+  RR = R; // Copies checkerboard for insert
+  
   typedef typename vobj::scalar_object sobj;
   typedef typename vobj::vector_type vector_type;
-
-  int Nblock = X.Grid()->GlobalDimensions()[Orthog];
-
-  GridBase *FullGrid  = X.Grid();
-  //  GridBase *SliceGrid = makeSubSliceGrid(FullGrid,Orthog);
-
-  //  Lattice<vobj> Xslice(SliceGrid);
-  //  Lattice<vobj> Rslice(SliceGrid);
-
-  assert( FullGrid->_simd_layout[Orthog]==1);
-  //  int nh =  FullGrid->_ndimension;
-  //  int nl = SliceGrid->_ndimension;
-  //  int nl = nh-1;
-
-  //FIXME package in a convenient iterator
-  //Should loop over a plane orthogonal to direction "Orthog"
-  int stride=FullGrid->_slice_stride[Orthog];
-  int block =FullGrid->_slice_block [Orthog];
-  int nblock=FullGrid->_slice_nblock[Orthog];
-  int ostride=FullGrid->_ostride[Orthog];
-
-  autoView( X_v, X, CpuRead);
-  autoView( Y_v, Y, CpuRead);
-  autoView( R_v, R, CpuWrite);
-  thread_region
-  {
-    Vector<vobj> s_x(Nblock);
-
-    thread_for_collapse_in_region(2, n,nblock, {
-     for(int b=0;b<block;b++){
-      int o  = n*stride + b;
-
-      for(int i=0;i<Nblock;i++){
-	s_x[i] = X_v[o+i*ostride];
-      }
-
-      vobj dot;
-      for(int i=0;i<Nblock;i++){
-	dot = Y_v[o+i*ostride];
-	for(int j=0;j<Nblock;j++){
-	  dot = dot + s_x[j]*(scale*aa(j,i));
-	}
-	R_v[o+i*ostride]=dot;
-      }
-    }});
+  int Nslice = X.Grid()->GlobalDimensions()[Orthog];
+  for(int i=0;i<Nslice;i++){
+    ExtractSlice(Ys,Y,i,Orthog);
+    ExtractSlice(Rs,R,i,Orthog);
+    Rs=Ys;
+    for(int j=0;j<Nslice;j++){
+      ExtractSlice(Xs,X,j,Orthog);
+      Rs = Rs + Xs*(scale*aa(j,i));
+    }
+    InsertSlice(Rs,RR,i,Orthog);
   }
+  R=RR; // Copy back handles arguments aliasing case
+  delete SliceGrid;
 };
 
 template<class vobj>
-static void sliceMulMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<vobj> &X,int Orthog,RealD scale=1.0) 
-{    
-  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::vector_type vector_type;
-
-  int Nblock = X.Grid()->GlobalDimensions()[Orthog];
-
-  GridBase *FullGrid  = X.Grid();
-  //  GridBase *SliceGrid = makeSubSliceGrid(FullGrid,Orthog);
-  //  Lattice<vobj> Xslice(SliceGrid);
-  //  Lattice<vobj> Rslice(SliceGrid);
-
-  assert( FullGrid->_simd_layout[Orthog]==1);
-  //  int nh =  FullGrid->_ndimension;
-  //  int nl = SliceGrid->_ndimension;
-  //  int nl=1;
-
-  //FIXME package in a convenient iterator
-  // thread_for2d_in_region
-  //Should loop over a plane orthogonal to direction "Orthog"
-  int stride=FullGrid->_slice_stride[Orthog];
-  int block =FullGrid->_slice_block [Orthog];
-  int nblock=FullGrid->_slice_nblock[Orthog];
-  int ostride=FullGrid->_ostride[Orthog];
-  autoView( R_v, R, CpuWrite);
-  autoView( X_v, X, CpuRead);
-  thread_region
-  {
-    std::vector<vobj> s_x(Nblock);
-
-
-    thread_for_collapse_in_region( 2 ,n,nblock,{
-    for(int b=0;b<block;b++){
-      int o  = n*stride + b;
-
-      for(int i=0;i<Nblock;i++){
-	s_x[i] = X_v[o+i*ostride];
-      }
-
-      vobj dot;
-      for(int i=0;i<Nblock;i++){
-	dot = s_x[0]*(scale*aa(0,i));
-	for(int j=1;j<Nblock;j++){
-	  dot = dot + s_x[j]*(scale*aa(j,i));
-	}
-	R_v[o+i*ostride]=dot;
-      }
-    }});
-  }
+static void sliceMulMatrix (Lattice<vobj> &R,Eigen::MatrixXcd &aa,const Lattice<vobj> &X,int Orthog,RealD scale=1.0)
+{
+  R=Zero();
+  sliceMaddMatrix(R,aa,X,R,Orthog,scale);
 };
 
 
 template<class vobj>
 static void sliceInnerProductMatrix(  Eigen::MatrixXcd &mat, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int Orthog) 
 {
+  GridBase *SliceGrid = makeSubSliceGrid(lhs.Grid(),Orthog);
+
+  Lattice<vobj> ls(SliceGrid);
+  Lattice<vobj> rs(SliceGrid);
+  
   typedef typename vobj::scalar_object sobj;
   typedef typename vobj::vector_type vector_type;
-  
-  GridBase *FullGrid  = lhs.Grid();
-  //  GridBase *SliceGrid = makeSubSliceGrid(FullGrid,Orthog);
-  
-  int Nblock = FullGrid->GlobalDimensions()[Orthog];
-  
-  //  Lattice<vobj> Lslice(SliceGrid);
-  //  Lattice<vobj> Rslice(SliceGrid);
-  
-  mat = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-
-  assert( FullGrid->_simd_layout[Orthog]==1);
-  //  int nh =  FullGrid->_ndimension;
-  //  int nl = SliceGrid->_ndimension;
-  //  int nl = nh-1;
-
-  //FIXME package in a convenient iterator
-  //Should loop over a plane orthogonal to direction "Orthog"
-  int stride=FullGrid->_slice_stride[Orthog];
-  int block =FullGrid->_slice_block [Orthog];
-  int nblock=FullGrid->_slice_nblock[Orthog];
-  int ostride=FullGrid->_ostride[Orthog];
-
-  typedef typename vobj::vector_typeD vector_typeD;
-
-  autoView( lhs_v, lhs, CpuRead);
-  autoView( rhs_v, rhs, CpuRead);
-  thread_region
-  {
-    std::vector<vobj> Left(Nblock);
-    std::vector<vobj> Right(Nblock);
-    Eigen::MatrixXcd  mat_thread = Eigen::MatrixXcd::Zero(Nblock,Nblock);
-
-    thread_for_collapse_in_region( 2, n,nblock,{
-    for(int b=0;b<block;b++){
-
-      int o  = n*stride + b;
-
-      for(int i=0;i<Nblock;i++){
-	Left [i] = lhs_v[o+i*ostride];
-	Right[i] = rhs_v[o+i*ostride];
-      }
-
-      for(int i=0;i<Nblock;i++){
-      for(int j=0;j<Nblock;j++){
-	auto tmp = innerProduct(Left[i],Right[j]);
-	auto rtmp = TensorRemove(tmp);
-	auto red  =  Reduce(rtmp);
-	mat_thread(i,j) += std::complex<double>(real(red),imag(red));
-      }}
-    }});
-    thread_critical
-    {
-      mat += mat_thread;
-    }  
+  int Nslice = lhs.Grid()->GlobalDimensions()[Orthog];
+  mat = Eigen::MatrixXcd::Zero(Nslice,Nslice);
+  for(int s=0;s<Nslice;s++){
+    ExtractSlice(ls,lhs,s,Orthog);
+    for(int ss=0;ss<Nslice;ss++){
+      ExtractSlice(rs,rhs,ss,Orthog);
+      mat(s,ss) = innerProduct(ls,rs);
+    }
   }
-
-  for(int i=0;i<Nblock;i++){
-  for(int j=0;j<Nblock;j++){
-    ComplexD sum = mat(i,j);
-    FullGrid->GlobalSum(sum);
-    mat(i,j)=sum;
-  }}
-
-  return;
+  delete SliceGrid;
 }
 
 NAMESPACE_END(Grid);

Grid/lattice/Lattice_transfer.h

@@ -981,8 +981,14 @@ void InsertSlice(const Lattice<vobj> &lowDim,Lattice<vobj> & higherDim,int slice
     hcoor[orthog] = slice;
     for(int d=0;d<nh;d++){
       if ( d!=orthog ) { 
-	hcoor[d]=lcoor[ddl++];
+	hcoor[d]=lcoor[ddl];
+	if ( hg->_checker_dim == d ) {
+	  hcoor[d]=hcoor[d]*2; // factor in the full coor for peekLocalSite
+	  lcoor[ddl]=lcoor[ddl]*2; // factor in the full coor for peekLocalSite
+	}
+	ddl++;
       }
+      
     }
     peekLocalSite(s,lowDimv,lcoor);
     pokeLocalSite(s,higherDimv,hcoor);
@@ -1003,6 +1009,7 @@ void ExtractSlice(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int slic
   assert(orthog<nh);
   assert(orthog>=0);
   assert(hg->_processors[orthog]==1);
+  lowDim.Checkerboard() = higherDim.Checkerboard();
 
   int dl; dl = 0;
   for(int d=0;d<nh;d++){
@@ -1020,11 +1027,16 @@ void ExtractSlice(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int slic
     Coordinate lcoor(nl);
     Coordinate hcoor(nh);
     lg->LocalIndexToLocalCoor(idx,lcoor);
-    int ddl=0;
     hcoor[orthog] = slice;
+    int ddl=0;
     for(int d=0;d<nh;d++){
       if ( d!=orthog ) { 
-	hcoor[d]=lcoor[ddl++];
+	hcoor[d]=lcoor[ddl];
+	if ( hg->_checker_dim == d ) {
+	  hcoor[d]=hcoor[d]*2;     // factor in the full gridd coor for peekLocalSite
+	  lcoor[ddl]=lcoor[ddl]*2; // factor in the full coor for peekLocalSite
+	}
+	ddl++;
       }
     }
     peekLocalSite(s,higherDimv,hcoor);

systems/Aurora/benchmarks/bench1.pbs

@@ -33,11 +33,11 @@ export MPICH_OFI_NIC_POLICY=GPU
 
 # 12 ppn, 2 nodes, 24 ranks
 #
-CMD="mpiexec -np 12 -ppn 12  -envall \
+CMD="mpiexec -np 1 -ppn 1  -envall \
 	     ./gpu_tile_compact.sh \
-	     ./Benchmark_comms_host_device --mpi 2.2.1.3 --grid 24.32.32.24 \
+	     ./Benchmark_usqcd --mpi 1.1.1.1 --grid 24.32.32.24 \
 		--shm-mpi 0 --shm 2048 --device-mem 32000 --accelerator-threads 32" 
-#$CMD | tee 1node.comms
+$CMD | tee usqcd.log
 
 
 CMD="mpiexec -np 1 -ppn 1  -envall \
@@ -50,7 +50,7 @@ CMD="mpiexec -np 12 -ppn 12  -envall \
 	     ./gpu_tile_compact.sh \
 	     ./Benchmark_dwf_fp32 --mpi 2.2.1.3 --grid 32.32.32.48 \
 		--shm-mpi 0 --shm 2048 --device-mem 32000 --accelerator-threads 32 --comms-overlap"
-$CMD | tee 1node.32.32.32.48.dwf
+#$CMD | tee 1node.32.32.32.48.dwf
 
 
 CMD="mpiexec -np 12 -ppn 12  -envall \

systems/Aurora/config-command

@@ -1,6 +1,6 @@
 
-export LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen -Xs -device -Xs pvc -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel  -fsycl  -lsycl " 
-export CXXFLAGS="-O3 -fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -qmkl=parallel  -fsycl -fno-exceptions -fsycl-targets=spir64_gen -Xs -device -Xs pvc "
+export LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel  -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel  -fsycl  -lsycl " 
+export CXXFLAGS="-O3 -fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -qmkl=parallel  -fsycl -fno-exceptions "
 ../../configure \
 	--enable-simd=GPU \
 	--enable-gen-simd-width=64 \

tests/debug/Test_general_coarse_hdcg.cc

@@ -77,13 +77,13 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  const int Ls=24;
-  const int nbasis = 60;
+  const int Ls=16;
+  const int nbasis = 40;
   const int cb = 0 ;
-  RealD mass=0.00078;
+  RealD mass=0.01;
   RealD M5=1.8;
-  RealD b=1.5;
-  RealD c=0.5;
+  RealD b=1.0;
+  RealD c=0.0;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
 								   GridDefaultSimd(Nd,vComplex::Nsimd()),
@@ -117,7 +117,7 @@ int main (int argc, char ** argv)
   LatticeGaugeField Umu(UGrid);
 
   FieldMetaData header;
-  std::string file("ckpoint_EODWF_lat.125");
+  std::string file("ckpoint_lat.4000");
   NerscIO::readConfiguration(Umu,header,file);
 
   //////////////////////// Fermion action //////////////////////////////////
@@ -155,7 +155,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
   std::cout << "Refine Subspace"<<std::endl;
   std::cout << "**************************************"<<std::endl;
-  Aggregates.RefineSubspace(HermOpEO,0.001,1.0e-3,3000); // 172 iters
+  //  Aggregates.RefineSubspace(HermOpEO,0.01,1.0e-3,1000); 
   
   std::cout << "**************************************"<<std::endl;
   std::cout << "Coarsen after refine"<<std::endl;
@@ -167,7 +167,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
   ConjugateGradient<CoarseVector>  coarseCG(4.0e-2,20000,true);
     
-  const int nrhs=12;
+  const int nrhs=8;
     
   Coordinate mpi=GridDefaultMpi();
   Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
@@ -185,7 +185,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
 
   typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
-  Chebyshev<CoarseVector>      IRLCheby(0.01,42.0,301);  // 1 iter
+  Chebyshev<CoarseVector>      IRLCheby(0.05,40.0,101);  // 1 iter
   MrhsHermMatrix MrhsCoarseOp     (mrhs);
 
   CoarseVector pm_src(CoarseMrhs);
@@ -193,8 +193,10 @@ int main (int argc, char ** argv)
   PowerMethod<CoarseVector>       cPM;
   cPM(MrhsCoarseOp,pm_src);
 
-  int Nk=192;
-  int Nm=384;
+  int Nk=nrhs;
+  int Nm=Nk*3;
+  int Nk=36;
+  int Nm=144;
   int Nstop=Nk;
   int Nconv_test_interval=1;
   
@@ -208,7 +210,7 @@ int main (int argc, char ** argv)
 							  nrhs,
 							  Nk,
 							  Nm,
-							  1e-5,10);
+							  1e-4,10);
 
   int Nconv;
   std::vector<RealD>            eval(Nm);
@@ -250,9 +252,9 @@ int main (int argc, char ** argv)
   // Extra HDCG parameters
   //////////////////////////
   int maxit=3000;
-  ConjugateGradient<CoarseVector>  CG(5.0e-2,maxit,false);
+  ConjugateGradient<CoarseVector>  CG(2.0e-1,maxit,false);
   RealD lo=2.0;
-  int ord = 7;
+  int ord = 9;
 
   DoNothingGuesser<CoarseVector> DoNothing;
   HPDSolver<CoarseVector> HPDSolveMrhs(MrhsCoarseOp,CG,DoNothing);

tests/debug/Test_general_coarse_hdcg_phys48.cc

@@ -287,8 +287,8 @@ int main (int argc, char ** argv)
   } else {
     //    Aggregates.CreateSubspaceMultishift(RNG5,HermOpEO,
     //					0.0003,1.0e-5,2000); // Lo, tol, maxit
-    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500); <== last run
-    //    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); // 176 with refinement
+    //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500); <== last run
+    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); // 176 with refinement
     //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.001,3000,1500,200,0.0); // Attempt to resurrect
     SaveBasis(Aggregates,subspace_file);
   }

tests/debug/Test_general_coarse_hdcg_phys48_lanczos.cc

@@ -0,0 +1,761 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_general_coarse_hdcg.cc
+
+    Copyright (C) 2023
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
+#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
+#include <Grid/algorithms/iterative/AdefMrhs.h>
+#include <Grid/algorithms/iterative/PowerSpectrum.h>
+#include <Grid/algorithms/iterative/BlockConjugateGradient.h>
+
+using namespace std;
+using namespace Grid;
+
+template<class aggregation>
+void SaveFineEvecs(aggregation &Agg,std::string file)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacWriter WR(Agg[0].Grid()->IsBoss());
+  WR.open(file);
+  for(int b=0;b<Agg.size();b++){
+    WR.writeScidacFieldRecord(Agg[b],record,0,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+  }
+  WR.close();
+#endif
+}
+template<class aggregation>
+void SaveBasis(aggregation &Agg,std::string file)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacWriter WR(Agg.FineGrid->IsBoss());
+  WR.open(file);
+  for(int b=0;b<Agg.subspace.size();b++){
+    WR.writeScidacFieldRecord(Agg.subspace[b],record,0,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+    //    WR.writeScidacFieldRecord(Agg.subspace[b],record);
+  }
+  WR.close();
+#endif
+}
+template<class aggregation>
+void LoadBasis(aggregation &Agg, std::string file)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacReader RD ;
+  RD.open(file);
+  for(int b=0;b<Agg.subspace.size();b++){
+    RD.readScidacFieldRecord(Agg.subspace[b],record,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+    //    RD.readScidacFieldRecord(Agg.subspace[b],record,0);
+  }    
+  RD.close();
+#endif
+}
+
+template<class aggregation>
+void LoadBasisSkip(aggregation &Agg, std::string file,int N,LatticeFermionF & tmp)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacReader RD ;
+  
+  RD.open(file);
+  for(int b=0;b<Agg.subspace.size();b++){
+    for(int n=0;n<N;n++){
+      RD.readScidacFieldRecord(tmp,record,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+      if(n==0) precisionChange(Agg.subspace[b],tmp);
+    }
+    //    RD.readScidacFieldRecord(Agg.subspace[b],record,0);
+  }    
+  RD.close();
+#endif
+}
+template<class aggregation>
+void LoadBasisSum(aggregation &Agg, std::string file,int N,LatticeFermionF & tmp)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacReader RD ;
+  
+  LatticeFermionF sum(tmp.Grid());
+  RD.open(file);
+  for(int b=0;b<Agg.subspace.size();b++){
+    sum=Zero();
+    for(int n=0;n<N;n++){
+      RD.readScidacFieldRecord(tmp,record,Grid::BinaryIO::BINARYIO_LEXICOGRAPHIC);
+      sum=sum+tmp;
+    }
+    precisionChange(Agg.subspace[b],sum);
+    //    RD.readScidacFieldRecord(Agg.subspace[b],record,0);
+  }    
+  RD.close();
+#endif
+}
+
+template<class CoarseVector>
+void SaveEigenvectors(std::vector<RealD>            &eval,
+		      std::vector<CoarseVector>     &evec,
+		      std::string evec_file,
+		      std::string eval_file)
+{
+#ifdef HAVE_LIME
+  emptyUserRecord record;
+  ScidacWriter WR(evec[0].Grid()->IsBoss());
+  WR.open(evec_file);
+  for(int b=0;b<evec.size();b++){
+    WR.writeScidacFieldRecord(evec[b],record,0,0);
+  }
+  WR.close();
+  XmlWriter WRx(eval_file);
+  write(WRx,"evals",eval);
+#endif
+}
+template<class CoarseVector>
+void LoadEigenvectors(std::vector<RealD>            &eval,
+		      std::vector<CoarseVector>     &evec,
+		      std::string evec_file,
+		      std::string eval_file)
+{
+#ifdef HAVE_LIME
+    XmlReader RDx(eval_file);
+    read(RDx,"evals",eval);
+    emptyUserRecord record;
+
+    Grid::ScidacReader RD ;
+    RD.open(evec_file);
+    assert(evec.size()==eval.size());
+    for(int k=0;k<eval.size();k++) {
+      RD.readScidacFieldRecord(evec[k],record);
+    }
+    RD.close();
+#endif
+}
+
+// Want Op in CoarsenOp to call MatPcDagMatPc
+template<class Field>
+class HermOpAdaptor : public LinearOperatorBase<Field>
+{
+  LinearOperatorBase<Field> & wrapped;
+public:
+  HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme)  {};
+  void Op     (const Field &in, Field &out)   { wrapped.HermOp(in,out);  }
+  void HermOp(const Field &in, Field &out)    { wrapped.HermOp(in,out); }
+  void AdjOp     (const Field &in, Field &out){ wrapped.HermOp(in,out);  }
+  void OpDiag (const Field &in, Field &out)                  {    assert(0);  }
+  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
+  void OpDirAll  (const Field &in, std::vector<Field> &out)  {    assert(0);  };
+  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
+};
+
+template<class Field> class FixedCGPolynomial : public LinearFunction<Field>
+{
+public:
+  using LinearFunction<Field>::operator();
+  typedef LinearOperatorBase<Field> FineOperator;
+  FineOperator   & _SmootherOperator;
+  ConjugateGradientPolynomial<Field>  CG;
+  int iters;
+  bool record;
+  int replay_count;
+  FixedCGPolynomial(int _iters, FineOperator &SmootherOperator) :
+    _SmootherOperator(SmootherOperator),
+    iters(_iters),
+    record(true),
+    CG(0.0,_iters,false)
+  {
+    std::cout << GridLogMessage<<" FixedCGPolynomial order "<<iters<<std::endl;
+    replay_count = 0;
+  };
+  void operator() (const Field &in, Field &out) 
+  {
+#if 1
+    GridBase *grid = in.Grid();
+    Field Mx0(grid);
+    Field r0(grid);
+    Field Minvr0(grid);
+
+    _SmootherOperator.HermOp(out,Mx0);
+
+    r0 = in - Mx0;
+
+    Minvr0 = Zero();
+    Minvr0.Checkerboard()=in.Checkerboard();
+    
+    if ( record ) {
+      std::cout << " FixedCGPolynomial recording polynomial "<<std::endl;
+      CG.Solve(_SmootherOperator,r0,Minvr0);
+      record = false;
+      /*
+      std::cout << "P(x) = 0 "<<std::endl;
+      for(int i=0;i<CG.polynomial.size();i++){
+	std::cout<<" + "<< CG.polynomial[i]<<" * (x**"<<i<<")"<<std::endl;
+	}
+      */
+      Field tmp(Minvr0.Grid());
+      CG.CGsequenceHermOp(_SmootherOperator,r0,tmp);
+      tmp = tmp - Minvr0;
+      std::cout << " CGsequence error "<<norm2(tmp)<<" / "<<norm2(out)<<std::endl;
+    } else {
+      std::cout << " FixedCGPolynomial replaying polynomial "<<std::endl;
+      CG.CGsequenceHermOp(_SmootherOperator,r0,Minvr0);
+      if ( replay_count %5== 0 ) record=true;
+      replay_count++;
+    }
+    out = out + Minvr0;
+    _SmootherOperator.HermOp(out,r0);
+    r0 = r0 - in;
+    RealD rr=norm2(r0);
+    RealD ss=norm2(in);
+    std::cout << " FixedCGPolynomial replayed polynomial resid "<<::sqrt(rr/ss)<<std::endl;
+#else
+    out = Zero();
+    out.Checkerboard()=in.Checkerboard();
+    if ( record ) {
+      std::cout << " FixedCGPolynomial recording polynomial "<<std::endl;
+      CG.Solve(_SmootherOperator,in,out);
+      record = false;
+      std::cout << "P(x) = 0 "<<std::endl;
+      for(int i=0;i<CG.polynomial.size();i++){
+	std::cout<<" + "<< CG.polynomial[i]<<" * (x**"<<i<<")"<<std::endl;
+      }
+      Field tmp(in.Grid());
+      CG.CGsequenceHermOp(_SmootherOperator,in,tmp);
+      tmp = tmp - out;
+      std::cout << " CGsequence error "<<norm2(tmp)<<" / "<<norm2(out)<<std::endl;
+    } else {
+      std::cout << " FixedCGPolynomial replaying polynomial "<<std::endl;
+      CG.CGsequenceHermOp(_SmootherOperator,in,out);
+      if ( replay_count %5== 5 ) record=true;
+      replay_count++;
+    }
+#endif
+    
+  }
+  void operator() (const std::vector<Field> &in, std::vector<Field> &out)
+  {
+    for(int i=0;i<out.size();i++){
+      out[i]=Zero();
+    }
+    int blockDim = 0;//not used for BlockCGVec
+    BlockConjugateGradient<Field>    BCGV  (BlockCGrQVec,blockDim,0.0,iters,false);
+    BCGV(_SmootherOperator,in,out);
+  }
+  
+};
+template<class Field> class CGSmoother : public LinearFunction<Field>
+{
+public:
+  using LinearFunction<Field>::operator();
+  typedef LinearOperatorBase<Field> FineOperator;
+  FineOperator   & _SmootherOperator;
+  int iters;
+  CGSmoother(int _iters, FineOperator &SmootherOperator) :
+    _SmootherOperator(SmootherOperator),
+    iters(_iters)
+  {
+    std::cout << GridLogMessage<<" Mirs smoother order "<<iters<<std::endl;
+  };
+  void operator() (const Field &in, Field &out) 
+  {
+    ConjugateGradient<Field>  CG(0.0,iters,false); // non-converge is just fine in a smoother
+
+    out=Zero();
+
+    CG(_SmootherOperator,in,out);
+  }
+};
+
+
+RealD InverseApproximation(RealD x){
+  return 1.0/x;
+}
+template<class Field> class ChebyshevSmoother : public LinearFunction<Field>
+{
+public:
+  using LinearFunction<Field>::operator();
+  typedef LinearOperatorBase<Field> FineOperator;
+  FineOperator   & _SmootherOperator;
+  Chebyshev<Field> Cheby;
+  ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
+    _SmootherOperator(SmootherOperator),
+    Cheby(_lo,_hi,_ord,InverseApproximation)
+  {
+    std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
+  };
+  void operator() (const Field &in, Field &out) 
+  {
+    //    Field r(out.Grid());
+    Cheby(_SmootherOperator,in,out);
+    //    _SmootherOperator.HermOp(out,r);
+    //    r=r-in;
+    //    RealD rr=norm2(r);
+    //    RealD ss=norm2(in);
+    //    std::cout << GridLogMessage<<" Chebyshev smoother resid "<<::sqrt(rr/ss)<<std::endl;
+  }
+};
+
+template<class Field> class ChebyshevInverter : public LinearFunction<Field>
+{
+public:
+  using LinearFunction<Field>::operator();
+  typedef LinearOperatorBase<Field> FineOperator;
+  FineOperator   & _Operator;
+  Chebyshev<Field> Cheby;
+  ChebyshevInverter(RealD _lo,RealD _hi,int _ord, FineOperator &Operator) :
+    _Operator(Operator),
+    Cheby(_lo,_hi,_ord,InverseApproximation)
+  {
+    std::cout << GridLogMessage<<" Chebyshev Inverter order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
+  };
+  void operator() (const Field &in, Field &out) 
+  {
+    Field r(in.Grid());
+    Field AinvR(in.Grid());
+    _Operator.HermOp(out,r);
+    r = in - r; // b - A x
+    Cheby(_Operator,r,AinvR); // A^{-1} ( b - A x ) ~ A^{-1} b - x
+    out = out + AinvR;
+    _Operator.HermOp(out,r);
+    r = in - r; // b - A x
+    RealD rr = norm2(r);
+    RealD ss = norm2(in);
+    std::cout << "ChebshevInverse resid " <<::sqrt(rr/ss)<<std::endl;
+  }
+};
+
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int sample=1;
+  if( GridCmdOptionExists(argv,argv+argc,"--sample") ){
+    std::string arg;
+    arg = GridCmdOptionPayload(argv,argv+argc,"--sample");
+    GridCmdOptionInt(arg,sample);
+  }
+  
+  const int Ls=24;
+  const int nbasis = 64;
+  const int cb = 0 ;
+  RealD mass=0.00078;
+
+  if( GridCmdOptionExists(argv,argv+argc,"--mass") ){
+    std::string arg;
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mass");
+    GridCmdOptionFloat(arg,mass);
+  }
+
+  RealD M5=1.8;
+  RealD b=1.5;
+  RealD c=0.5;
+
+  std::cout << GridLogMessage << " *************************** " <<std::endl;
+  std::cout << GridLogMessage << " Mass " <<mass<<std::endl;
+  std::cout << GridLogMessage << " M5   " <<M5<<std::endl;
+  std::cout << GridLogMessage << " Ls   " <<Ls<<std::endl;
+  std::cout << GridLogMessage << " b    " <<b<<std::endl;
+  std::cout << GridLogMessage << " c    " <<c<<std::endl;
+  std::cout << GridLogMessage << " *************************** " <<std::endl;
+  
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  //////////////////////////////////////////
+  // Single precision grids -- lanczos + smoother
+  //////////////////////////////////////////
+  GridCartesian         * UGridF   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
+								   GridDefaultSimd(Nd,vComplexF::Nsimd()),
+								   GridDefaultMpi());
+  GridRedBlackCartesian * UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF);
+  GridCartesian         * FGridF   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF);
+  GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
+
+  ///////////////////////// Configuration /////////////////////////////////
+  LatticeGaugeField Umu(UGrid);
+
+  FieldMetaData header;
+  std::string file("ckpoint_lat.1000");
+  NerscIO::readConfiguration(Umu,header,file);
+
+  //////////////////////// Fermion action //////////////////////////////////
+  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
+  SchurDiagMooeeOperator<MobiusFermionD, LatticeFermion> HermOpEO(Ddwf);
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "         Fine Power method            "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+
+  {
+    LatticeFermionD pm_src(FrbGrid);
+    pm_src = ComplexD(1.0);
+    PowerMethod<LatticeFermionD>       fPM;
+    fPM(HermOpEO,pm_src);
+  }
+
+  if(0)
+  {
+
+    std::cout << "**************************************"<<std::endl;
+    std::cout << "         Fine Lanczos           "<<std::endl;
+    std::cout << "**************************************"<<std::endl;
+
+    typedef LatticeFermionF FermionField;
+    LatticeGaugeFieldF UmuF(UGridF);
+    precisionChange(UmuF,Umu);
+    MobiusFermionF DdwfF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,b,c);
+    SchurDiagMooeeOperator<MobiusFermionF, LatticeFermionF> HermOpEOF(DdwfF);
+
+    const int Fine_Nstop = 200;
+    const int Fine_Nk = 200;
+    const int Fine_Np = 200;
+    const int Fine_Nm = Fine_Nk+Fine_Np;
+    const int Fine_MaxIt= 10;
+
+    RealD Fine_resid = 1.0e-4;
+    std::cout << GridLogMessage << "Fine Lanczos "<<std::endl;
+    std::cout << GridLogMessage << "Nstop "<<Fine_Nstop<<std::endl;
+    std::cout << GridLogMessage << "Nk "<<Fine_Nk<<std::endl;
+    std::cout << GridLogMessage << "Np "<<Fine_Np<<std::endl;
+    std::cout << GridLogMessage << "resid "<<Fine_resid<<std::endl;
+
+    Chebyshev<FermionField> Cheby(0.002,92.0,401);
+    //    Chebyshev<FermionField> Cheby(0.1,92.0,401);
+    FunctionHermOp<FermionField> OpCheby(Cheby,HermOpEOF);
+    PlainHermOp<FermionField> Op     (HermOpEOF);
+    ImplicitlyRestartedLanczos<FermionField> IRL(OpCheby,Op,Fine_Nstop,Fine_Nk,Fine_Nm,Fine_resid,Fine_MaxIt);
+    std::vector<RealD>          Fine_eval(Fine_Nm);
+    FermionField                Fine_src(FrbGridF); 
+    Fine_src = ComplexF(1.0);
+    std::vector<FermionField> Fine_evec(Fine_Nm,FrbGridF);
+
+    int Fine_Nconv;
+    std::cout << GridLogMessage <<" Calling IRL.calc single prec"<<std::endl;
+    IRL.calc(Fine_eval,Fine_evec,Fine_src,Fine_Nconv);
+
+    std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.evecF");
+    SaveFineEvecs(Fine_evec,evec_file);
+  }
+
+
+  //////////////////////////////////////////
+  // Construct a coarsened grid with 4^4 cell
+  //////////////////////////////////////////
+  Coordinate Block({4,4,6,4});
+  Coordinate clatt = GridDefaultLatt();
+  for(int d=0;d<clatt.size();d++){
+    clatt[d] = clatt[d]/Block[d];
+  }
+
+  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt,
+							    GridDefaultSimd(Nd,vComplex::Nsimd()),
+							    GridDefaultMpi());;
+  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+
+  ///////////////////////// RNGs /////////////////////////////////
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+  std::vector<int> cseeds({5,6,7,8});
+
+  GridParallelRNG          RNG5(FGrid);   RNG5.SeedFixedIntegers(seeds5);
+  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
+
+  
+  typedef HermOpAdaptor<LatticeFermionD> HermFineMatrix;
+  HermFineMatrix FineHermOp(HermOpEO);
+
+  ////////////////////////////////////////////////////////////
+  ///////////// Coarse basis and Little Dirac Operator ///////
+  ////////////////////////////////////////////////////////////
+  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
+  typedef LittleDiracOperator::CoarseVector CoarseVector;
+
+  NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
+
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
+  Subspace Aggregates(Coarse5d,FrbGrid,cb);
+
+  ////////////////////////////////////////////////////////////
+  // Need to check about red-black grid coarsening
+  ////////////////////////////////////////////////////////////
+  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.mixed.2500.60");
+  //  //  std::string subspace_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.new.62");
+  std::string refine_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Subspace.phys48.evecF");
+  //  std::string refine_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/Refine.phys48.mixed.2500.60");
+  std::string ldop_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/LittleDiracOp.phys48.mixed.60");
+  std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/evecs.scidac");
+  std::string eval_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/eval.xml");
+  bool load_agg=true;
+  bool load_refine=true;
+  bool load_mat=false;
+  bool load_evec=false;
+
+  int refine=1;
+  if ( load_agg ) {
+    if ( !(refine) || (!load_refine) ) { 
+      LoadBasis(Aggregates,subspace_file);
+    }
+  } else {
+    //    Aggregates.CreateSubspaceMultishift(RNG5,HermOpEO,
+    //    					0.0003,1.0e-5,2000); // Lo, tol, maxit
+    //    Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,95.,0.01,1500);// <== last run
+    Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); 
+    SaveBasis(Aggregates,subspace_file);
+  }
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "Building MultiRHS Coarse operator"<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  ConjugateGradient<CoarseVector>  coarseCG(4.0e-2,20000,true);
+    
+  const int nrhs=12;
+    
+  Coordinate mpi=GridDefaultMpi();
+  Coordinate rhMpi ({1,1,mpi[0],mpi[1],mpi[2],mpi[3]});
+  Coordinate rhLatt({nrhs,1,clatt[0],clatt[1],clatt[2],clatt[3]});
+  Coordinate rhSimd({vComplex::Nsimd(),1, 1,1,1,1});
+    
+  GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); 
+  typedef MultiGeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> MultiGeneralCoarsenedMatrix_t;
+  MultiGeneralCoarsenedMatrix_t mrhs(geom,CoarseMrhs);
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "         Coarse Lanczos               "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+
+  typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
+  Chebyshev<CoarseVector>      IRLCheby(0.005,42.0,301);  // 1 iter
+  MrhsHermMatrix MrhsCoarseOp     (mrhs);
+
+  //  CoarseVector pm_src(CoarseMrhs);
+  //  pm_src = ComplexD(1.0);
+  //  PowerMethod<CoarseVector>       cPM;   cPM(MrhsCoarseOp,pm_src);
+
+  int Nk=192;
+  int Nm=384;
+  int Nstop=Nk;
+  int Nconv_test_interval=1;
+  
+  ImplicitlyRestartedBlockLanczosCoarse<CoarseVector> IRL(MrhsCoarseOp,
+							  Coarse5d,
+							  CoarseMrhs,
+							  nrhs,
+							  IRLCheby,
+							  Nstop,
+							  Nconv_test_interval,
+							  nrhs,
+							  Nk,
+							  Nm,
+							  1e-5,10);
+
+  int Nconv;
+  std::vector<RealD>            eval(Nm);
+  std::vector<CoarseVector>     evec(Nm,Coarse5d);
+  std::vector<CoarseVector>     c_src(nrhs,Coarse5d);
+
+  ///////////////////////
+  // Deflation guesser object
+  ///////////////////////
+  MultiRHSDeflation<CoarseVector> MrhsGuesser;
+
+  //////////////////////////////////////////
+  // Block projector for coarse/fine
+  //////////////////////////////////////////
+  MultiRHSBlockProject<LatticeFermionD> MrhsProjector;
+
+  //////////////////////////
+  // Extra HDCG parameters
+  //////////////////////////
+  int maxit=300;
+  ConjugateGradient<CoarseVector>  CG(5.0e-2,maxit,false);
+  ConjugateGradient<CoarseVector>  CGstart(5.0e-2,maxit,false);
+  RealD lo=2.0;
+  int ord = 7;
+  //  int ord = 11;
+
+  int blockDim = 0;//not used for BlockCG
+  BlockConjugateGradient<CoarseVector>    BCG  (BlockCGrQ,blockDim,5.0e-5,maxit,true);
+
+  DoNothingGuesser<CoarseVector> DoNothing;
+  //  HPDSolver<CoarseVector> HPDSolveMrhs(MrhsCoarseOp,CG,DoNothing);
+  //  HPDSolver<CoarseVector> HPDSolveMrhsStart(MrhsCoarseOp,CGstart,DoNothing);
+  //  HPDSolver<CoarseVector> HPDSolveMrhs(MrhsCoarseOp,BCG,DoNothing);
+  //  HPDSolver<CoarseVector> HPDSolveMrhsRefine(MrhsCoarseOp,BCG,DoNothing);
+  //  FixedCGPolynomial<CoarseVector>  HPDSolveMrhs(maxit,MrhsCoarseOp);
+
+  ChebyshevInverter<CoarseVector> HPDSolveMrhs(1.0e-2,40.0,120,MrhsCoarseOp);  //
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(1.0e-2,40.0,110,MrhsCoarseOp);  // 114 iter with Chebysmooth and BlockCG
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(1.0e-2,40.0,120,MrhsCoarseOp); // 138 iter with Chebysmooth
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(1.0e-2,40.0,200,MrhsCoarseOp); // 139 iter
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(3.0e-3,40.0,200,MrhsCoarseOp); // 137 iter, CG smooth, flex
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(1.0e-3,40.0,200,MrhsCoarseOp); // 146 iter, CG smooth, flex
+  //  ChebyshevInverter<CoarseVector> HPDSolveMrhs(3.0e-4,40.0,200,MrhsCoarseOp); // 156 iter, CG smooth, flex
+
+  /////////////////////////////////////////////////
+  // Mirs smoother
+  /////////////////////////////////////////////////
+  ShiftedHermOpLinearOperator<LatticeFermionD> ShiftedFineHermOp(HermOpEO,lo);
+  //  FixedCGPolynomial<LatticeFermionD> CGsmooth(ord,ShiftedFineHermOp) ;
+  //  CGSmoother<LatticeFermionD> CGsmooth(ord,ShiftedFineHermOp) ;
+  ChebyshevSmoother<LatticeFermionD> CGsmooth(2.0,92.0,8,HermOpEO) ;
+  
+  if ( load_refine ) {
+    //LoadBasis(Aggregates,refine_file);
+    LatticeFermionF conv_tmp(FrbGridF);
+    LoadBasisSum(Aggregates,refine_file,sample,conv_tmp);
+  } else {
+    Aggregates.RefineSubspace(HermOpEO,0.001,1.0e-3,3000); // 172 iters
+    SaveBasis(Aggregates,refine_file);
+  }
+  Aggregates.Orthogonalise();
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "Coarsen after refine"<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  mrhs.CoarsenOperator(FineHermOp,Aggregates,Coarse5d);
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << " Recompute coarse evecs  "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  evec.resize(Nm,Coarse5d);
+  eval.resize(Nm);
+  for(int r=0;r<nrhs;r++){
+    random(CRNG,c_src[r]);
+  }
+ IRL.calc(eval,evec,c_src,Nconv,LanczosType::irbl);
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << " Reimport coarse evecs  "<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  MrhsGuesser.ImportEigenBasis(evec,eval);
+
+  std::cout << "**************************************"<<std::endl;
+  std::cout << " Setting up mRHS HDCG"<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  MrhsProjector.Allocate(nbasis,FrbGrid,Coarse5d);
+  MrhsProjector.ImportBasis(Aggregates.subspace);
+      
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "Calling mRHS HDCG"<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+  TwoLevelADEF2mrhs<LatticeFermion,CoarseVector>
+    HDCGmrhs(1.0e-8, 300,
+	     FineHermOp,
+	     CGsmooth,
+	     HPDSolveMrhs,    // Used in M1
+	     HPDSolveMrhs,          // Used in Vstart
+	     MrhsProjector,
+	     MrhsGuesser,
+	     CoarseMrhs);
+    
+  std::vector<LatticeFermionD> src_mrhs(nrhs,FrbGrid);
+  std::vector<LatticeFermionD> res_mrhs(nrhs,FrbGrid);
+  LatticeFermionD result_accurate(FrbGrid);
+  LatticeFermionD result_sloppy(FrbGrid);
+  LatticeFermionD error(FrbGrid);
+  LatticeFermionD residual(FrbGrid);
+
+  for(int r=0;r<nrhs;r++){
+    random(RNG5,src_mrhs[r]);
+    res_mrhs[r]=Zero();
+  }
+  HDCGmrhs(src_mrhs,res_mrhs);
+  result_accurate = res_mrhs[0];
+
+#if 0
+  std::vector<RealD> tols({1.0e-3,1.0e-4,1.0e-5});
+
+  std::vector<RealD>   bins({1.0e-3,1.0e-2,1.0e-1,1.0,10.0,100.0});
+  std::vector<int>   orders({6000  ,4000  ,1000  ,500,500 ,500});
+  PowerSpectrum GraphicEqualizer(bins,orders);
+
+  for(auto tol : tols) {
+    
+    TwoLevelADEF2mrhs<LatticeFermion,CoarseVector>
+      HDCGmrhsSloppy(tol, 500,
+		     FineHermOp,
+		     CGsmooth,
+		     HPDSolveMrhs,    // Used in M1
+		     HPDSolveMrhs,    // Used in Vstart
+		     MrhsProjector,
+		     MrhsGuesser,
+		     CoarseMrhs);
+  
+    //  Solve again to 10^-5
+    for(int r=0;r<nrhs;r++){
+      res_mrhs[r]=Zero();
+    }
+    HDCGmrhsSloppy(src_mrhs,res_mrhs);
+    
+    result_sloppy = res_mrhs[0];
+    error = result_sloppy - result_accurate;
+    FineHermOp.HermOp(result_sloppy,residual);
+    residual = residual - src_mrhs[0];
+    
+    std::cout << "**************************************"<<std::endl;
+    std::cout << GridLogMessage << " Converged to tolerance "<< tol<<std::endl;
+    std::cout << GridLogMessage << " Absolute error "<<norm2(error)<<std::endl;
+    std::cout << GridLogMessage << " Residual       "<<norm2(residual)<<std::endl;
+    std::cout << "**************************************"<<std::endl;
+
+    std::cout << "**************************************"<<std::endl;
+    std::cout << GridLogMessage << " PowerSpectrum of error   "<<std::endl;
+    std::cout << "**************************************"<<std::endl;
+    GraphicEqualizer(FineHermOp,error);
+    std::cout << "**************************************"<<std::endl;
+    std::cout << GridLogMessage << " PowerSpectrum of residual   "<<std::endl;
+    std::cout << "**************************************"<<std::endl;
+    GraphicEqualizer(FineHermOp,residual);
+
+  };
+#endif
+  
+  // Standard CG
+#if 0
+  {
+  std::cout << "**************************************"<<std::endl;
+  std::cout << "Calling red black CG"<<std::endl;
+  std::cout << "**************************************"<<std::endl;
+      
+    LatticeFermion result(FrbGrid); result=Zero();
+    LatticeFermion    src(FrbGrid); random(RNG5,src);
+    result=Zero();
+
+    ConjugateGradient<LatticeFermionD>  CGfine(1.0e-8,30000,false);
+    CGfine(HermOpEO, src, result);
+  }
+#endif  
+  Grid_finalize();
+  return 0;
+}