Merge 32e6d58356 into 41d8adca95

BLAS_benchmark/compile-command

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

Grid/algorithms/blas/BatchedBlas.h

@@ -208,9 +208,6 @@ public:
     assert(Bkn.size()==batchCount);
     assert(Cmn.size()==batchCount);
 
-    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
@@ -270,6 +267,7 @@ public:
     assert(err==CUBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_SYCL
+    std::cerr << " Calling SYCL batched ZGEMM "<<std::endl;
       int64_t m64=m;
       int64_t n64=n;
       int64_t k64=k;
@@ -277,20 +275,10 @@ public:
       int64_t ldb64=ldb;
       int64_t ldc64=ldc;
       int64_t batchCount64=batchCount;
-
+      oneapi::mkl::transpose notransp =oneapi::mkl::transpose::N;
-      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,
+						  &notransp,
-						  &iOpB,
+						  &notransp,
 						  &m64,&n64,&k64,
 						  (ComplexD *) &alpha_p[0],
 						  (const ComplexD **)&Amk[0], (const int64_t *)&lda64,
@@ -299,100 +287,42 @@ public:
 						  (ComplexD **)&Cmn[0], (const int64_t *)&ldc64,
 						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
       synchronise();
-#if 0
+    std::cerr << " Called SYCL batched ZGEMM "<<std::endl;
-      // This code was used to check the mat mul on Sunspot/OneMKL
-      std::cerr << " Called SYCL batched ZGEMM OpA "<< OpA << " OpB "<<OpB <<std::endl;
       std::vector<ComplexD> A(m*k);  // pointer list to matrices
       std::vector<ComplexD> B(k*n);
       std::vector<ComplexD> C(m*n);
-      //      int sda = lda*k;
+      int sda = lda*k;
-      //      int sdb = ldb*k;
+      int sdb = ldb*k;
-      //      int sdc = ldc*n;
+      int sdc = ldc*n;
-      std::cerr << " Checking the GEMM results "<<std::endl;
       for (int p = 0; p < 1; ++p) {
-	ComplexD * Amk_p;  // pointer list to matrices
+	acceleratorCopyFromDevice((void *)&Amk[p][0],(void *)&A[0],m*k*sizeof(ComplexD));
-	ComplexD * Bkn_p;  // pointer list to matrices
+	acceleratorCopyFromDevice((void *)&Bkn[p][0],(void *)&B[0],k*n*sizeof(ComplexD));
-	ComplexD * Cmn_p;  // pointer list to matrices
+	acceleratorCopyFromDevice((void *)&Cmn[p][0],(void *)&C[0],m*n*sizeof(ComplexD));
-	acceleratorCopyFromDevice((void *)&Amk[p],(void *)&Amk_p,sizeof(ComplexD*));
-	acceleratorCopyFromDevice((void *)&Bkn[p],(void *)&Bkn_p,sizeof(ComplexD*));
-	acceleratorCopyFromDevice((void *)&Cmn[p],(void *)&Cmn_p,sizeof(ComplexD*));
-	std::cerr << " p " << p << " copied pointers "<<std::endl;
-	acceleratorCopyFromDevice((void *)Amk_p,(void *)&A[0],m*k*sizeof(ComplexD));
-	acceleratorCopyFromDevice((void *)Bkn_p,(void *)&B[0],k*n*sizeof(ComplexD));
-	acceleratorCopyFromDevice((void *)Cmn_p,(void *)&C[0],m*n*sizeof(ComplexD));
-	std::cerr << " p " << p << " copied matrices "<<std::endl;
-	std::cerr << " C[0] "<<C[0]<<std::endl;
-	std::cerr << " A[0] "<<A[0]<<std::endl;
-	std::cerr << " B[0] "<<B[0]<<std::endl;
-	std::cerr << " m "<<m<<std::endl;
-	std::cerr << " n "<<n<<std::endl;
-	std::cerr << " k "<<k<<std::endl;
 	for (int mm = 0; mm < m; ++mm) {
 	  for (int nn = 0; nn < n; ++nn) {
 	    ComplexD c_mn(0.0);
-	    for (int kk = 0; kk < k; ++kk) {
+	    for (int kk = 0; kk < k; ++kk)
-	      int idx_a, idx_b;
+	      c_mn += A[mm + kk*lda ] * B[kk + nn*ldb];
-	      //    int lda = m; // m x k column major
+	    std::cout << " beta "<<beta<<" C_"<<mm<<","<<nn<<" "<<c_mn<<" "<<C[mm + nn*ldc]<<std::endl;
-	      //    int ldb = k; // k x n column major
-	      //    int ldc = m; // m x b column major
-	      if(OpA!=GridBLAS_OP_N) {
-		idx_a =kk + mm*lda;
-	      } else {
-		idx_a =mm + kk*lda;
-	      }
-	      if(OpB!=GridBLAS_OP_N) {
-		idx_b =nn + kk*ldb;
-	      } else {
-		idx_b =kk + nn*ldb;
-	      }
-	      //	      std::cerr << " idx_a "<<idx_a<<" idx_b "<<idx_b<<std::endl;
-
-	      ComplexD Ac = A[idx_a];
-	      ComplexD Bc = B[idx_b];
-	      if(OpA==GridBLAS_OP_C) Ac = conjugate(Ac);
-	      if(OpB==GridBLAS_OP_C) Bc = conjugate(Bc);
-	      
-	      c_mn += Ac*Bc;
-	    }
-	    std::cerr << " beta "<<beta<<" alpha "<<alpha<<" C_"<<mm<<","<<nn<<" "<<c_mn<<" "<<C[mm + nn*ldc]<<std::endl;
 	  }
 	}
       }
 #endif
-#endif
 #if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
+    // Need a default/reference implementation
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+    int sda = lda*k;
-	thread_for (p, batchCount, {
+    int sdb = ldb*k;
-	  Eigen::Map<Eigen::MatrixXcd> eAmk(Amk[p],m,k);
+    int sdc = ldc*n;
-	  Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn[p],k,n);
+    for (int p = 0; p < batchCount; ++p) {
-	  Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn[p],m,n);
+      for (int mm = 0; mm < m; ++mm) {
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+	for (int nn = 0; nn < n; ++nn) {
-        });
+	  ComplexD c_mn(0.0);
-      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+	  for (int kk = 0; kk < k; ++kk)
-	thread_for (p, batchCount, {
+	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Eigen::Map<Eigen::MatrixXcd> eAmk(Amk[p],k,m);
+	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	  Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn[p],k,n);
+	}
-	  Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXcd> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint() ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXcd> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXcd> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXcd> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint() ;
-	  } );
-      } else { 
-	assert(0);
       }
+    }
 #endif
      RealD t1=usecond();
      RealD flops = 8.0*m*n*k*batchCount;
@@ -414,9 +344,6 @@ public:
     RealD t2=usecond();
     int32_t batchCount = Amk.size();
 
-    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
@@ -484,20 +411,10 @@ public:
       int64_t ldb64=ldb;
       int64_t ldc64=ldc;
       int64_t batchCount64=batchCount;
-
+      oneapi::mkl::transpose notransp =oneapi::mkl::transpose::N;
-      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,
+						  &notransp,
-						  &iOpB,
+						  &notransp,
 						  &m64,&n64,&k64,
 						  (ComplexF *) &alpha_p[0],
 						  (const ComplexF **)&Amk[0], (const int64_t *)&lda64,
@@ -508,38 +425,22 @@ public:
     synchronise();
 #endif
 #if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
+    int sda = lda*k;
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+    int sdb = ldb*k;
-	thread_for (p, batchCount, {
+    int sdc = ldc*n;
-	  Eigen::Map<Eigen::MatrixXcf> eAmk(Amk[p],m,k);
+    ComplexF alphaf(real(alpha),imag(alpha));
-	  Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn[p],k,n);
+    ComplexF betaf(real(beta),imag(beta));
-	  Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn[p],m,n);
+    // Need a default/reference implementation
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+    for (int p = 0; p < batchCount; ++p) {
-	  });
+      for (int mm = 0; mm < m; ++mm) {
-      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_N) ) {
+	for (int nn = 0; nn < n; ++nn) {
-	thread_for (p, batchCount, {
+	  ComplexF c_mn(0.0);
-	  Eigen::Map<Eigen::MatrixXcf> eAmk(Amk[p],k,m);
+	  for (int kk = 0; kk < k; ++kk)
-	  Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn[p],k,n);
+	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn[p],m,n);
+	  Cmn[p][mm + nn*ldc] =  (alphaf)*c_mn + (betaf)*Cmn[p][mm + nn*ldc ];
-	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn ;
+	}
-	  });
-      } else if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_C) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXcf> eAmk(Amk[p],m,k);
-	  Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn.adjoint() ;
-	  });
-      } else if ( (OpA == GridBLAS_OP_C ) && (OpB == GridBLAS_OP_C) ) {
-	thread_for (p, batchCount, {
-	  Eigen::Map<Eigen::MatrixXcf> eAmk(Amk[p],k,m);
-	  Eigen::Map<Eigen::MatrixXcf> eBkn(Bkn[p],n,k);
-	  Eigen::Map<Eigen::MatrixXcf> eCmn(Cmn[p],m,n);
-	  eCmn = beta * eCmn + alpha * eAmk.adjoint() * eBkn.adjoint() ;
-	  } );
-      } else { 
-	assert(0);
       }
+    }
 #endif
      RealD t1=usecond();
      RealD flops = 8.0*m*n*k*batchCount;
@@ -562,9 +463,6 @@ public:
     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
@@ -631,20 +529,10 @@ public:
       int64_t ldb64=ldb;
       int64_t ldc64=ldc;
       int64_t batchCount64=batchCount;
-
+      oneapi::mkl::transpose notransp =oneapi::mkl::transpose::N;
-      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,
+						  &notransp,
-						  &iOpB,
+						  &notransp,
 						  &m64,&n64,&k64,
 						  (float *) &alpha_p[0],
 						  (const float **)&Amk[0], (const int64_t *)&lda64,
@@ -652,41 +540,23 @@ public:
 						  (float *) &beta_p[0],
 						  (float **)&Cmn[0], (const int64_t *)&ldc64,
 						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
-      synchronise();
+    synchronise();
 #endif
 #if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
+    int sda = lda*k;
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+    int sdb = ldb*k;
-	thread_for (p, batchCount, {
+    int sdc = ldc*n;
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],m,k);
+    // Need a default/reference implementation
-	  Eigen::Map<Eigen::MatrixXf> eBkn(Bkn[p],k,n);
+    for (int p = 0; p < batchCount; ++p) {
-	  Eigen::Map<Eigen::MatrixXf> eCmn(Cmn[p],m,n);
+      for (int mm = 0; mm < m; ++mm) {
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+	for (int nn = 0; nn < n; ++nn) {
-	  });
+	  RealD c_mn(0.0);
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
+	  for (int kk = 0; kk < k; ++kk)
-	thread_for (p, batchCount, {
+	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Eigen::Map<Eigen::MatrixXf> eAmk(Amk[p],k,m);
+	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	  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;
@@ -697,6 +567,7 @@ public:
   ///////////////////////////////////////////////////////////////////////////
   // Double precision real GEMM
   ///////////////////////////////////////////////////////////////////////////
+
   void gemmBatched(GridBLASOperation_t OpA,
 		   GridBLASOperation_t OpB,
 		   int m,int n, int k,
@@ -709,9 +580,6 @@ public:
     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
@@ -779,20 +647,10 @@ public:
       int64_t ldb64=ldb;
       int64_t ldc64=ldc;
       int64_t batchCount64=batchCount;
-
+      oneapi::mkl::transpose notransp =oneapi::mkl::transpose::N;
-      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,
+						  &notransp,
-						  &iOpB,
+						  &notransp,
 						  &m64,&n64,&k64,
 						  (double *) &alpha_p[0],
 						  (const double **)&Amk[0], (const int64_t *)&lda64,
@@ -800,96 +658,144 @@ public:
 						  (double *) &beta_p[0],
 						  (double **)&Cmn[0], (const int64_t *)&ldc64,
 						  (int64_t)1,&batchCount64,std::vector<sycl::event>());
-      synchronise();
+    synchronise();
 #endif
 #if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation; use Eigen
+    int sda = lda*k;
-      if ( (OpA == GridBLAS_OP_N ) && (OpB == GridBLAS_OP_N) ) {
+    int sdb = ldb*k;
-	thread_for (p, batchCount, {
+    int sdc = ldc*n;
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],m,k);
+    // Need a default/reference implementation
-	  Eigen::Map<Eigen::MatrixXd> eBkn(Bkn[p],k,n);
+    for (int p = 0; p < batchCount; ++p) {
-	  Eigen::Map<Eigen::MatrixXd> eCmn(Cmn[p],m,n);
+      for (int mm = 0; mm < m; ++mm) {
-	  eCmn = beta * eCmn + alpha * eAmk * eBkn ;
+	for (int nn = 0; nn < n; ++nn) {
-	  });
+	  RealD c_mn(0.0);
-      } else if ( (OpA == GridBLAS_OP_T ) && (OpB == GridBLAS_OP_N) ) {
+	  for (int kk = 0; kk < k; ++kk)
-	thread_for (p, batchCount, {
+	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Eigen::Map<Eigen::MatrixXd> eAmk(Amk[p],k,m);
+	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	  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;
   }
+  
+
+  
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  // Strided case used by benchmark, but generally unused in Grid
+  // Keep a code example in double complex, but don't generate the single and real variants for now
+  ////////////////////////////////////////////////////////////////////////////////////////////////
+  
+  void gemmStridedBatched(int m,int n, int k,
+			  ComplexD alpha,
+			  ComplexD* Amk,  // pointer list to matrices
+			  ComplexD* Bkn,
+			  ComplexD beta,
+			  ComplexD* Cmn,
+			  int batchCount)
+  {
+    // Use C-row major storage, so transpose calls
+    int lda = m; // m x k column major
+    int ldb = k; // k x n column major
+    int ldc = m; // m x b column major
+    int sda = m*k;
+    int sdb = k*n;
+    int sdc = m*n;
+    deviceVector<ComplexD> alpha_p(1);
+    deviceVector<ComplexD> beta_p(1);
+    acceleratorCopyToDevice((void *)&alpha,(void *)&alpha_p[0],sizeof(ComplexD));
+    acceleratorCopyToDevice((void *)&beta ,(void *)&beta_p[0],sizeof(ComplexD));
+    //    std::cout << "blasZgemmStridedBatched mnk  "<<m<<","<<n<<","<<k<<" count "<<batchCount<<std::endl;
+    //    std::cout << "blasZgemmStridedBatched ld   "<<lda<<","<<ldb<<","<<ldc<<std::endl;
+    //    std::cout << "blasZgemmStridedBatched sd   "<<sda<<","<<sdb<<","<<sdc<<std::endl;
+#ifdef GRID_HIP
+    auto err = hipblasZgemmStridedBatched(gridblasHandle,
+					  HIPBLAS_OP_N,
+					  HIPBLAS_OP_N,
+					  m,n,k,
+					  (hipblasDoubleComplex *) &alpha_p[0],
+					  (hipblasDoubleComplex *) Amk, lda, sda,
+					  (hipblasDoubleComplex *) Bkn, ldb, sdb,
+					  (hipblasDoubleComplex *) &beta_p[0],
+					  (hipblasDoubleComplex *) Cmn, ldc, sdc,
+					  batchCount);
+    assert(err==HIPBLAS_STATUS_SUCCESS);
+#endif
+#ifdef GRID_CUDA
+    cublasZgemmStridedBatched(gridblasHandle,
+			      CUBLAS_OP_N,
+			      CUBLAS_OP_N,
+			      m,n,k,
+			      (cuDoubleComplex *) &alpha_p[0],
+			      (cuDoubleComplex *) Amk, lda, sda,
+			      (cuDoubleComplex *) Bkn, ldb, sdb,
+			      (cuDoubleComplex *) &beta_p[0],
+			      (cuDoubleComplex *) Cmn, ldc, sdc,
+			      batchCount);
+#endif
+#if defined(GRID_SYCL) || defined(GRID_ONE_MKL)
+    oneapi::mkl::blas::column_major::gemm_batch(*gridblasHandle,
+						oneapi::mkl::transpose::N,
+						oneapi::mkl::transpose::N,
+						m,n,k,
+						alpha,
+						(const ComplexD *)Amk,lda,sda,
+						(const ComplexD *)Bkn,ldb,sdb,
+						beta,
+						(ComplexD *)Cmn,ldc,sdc,
+						batchCount);
+    synchronise();
+#endif
+#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP) && !defined(GRID_ONE_MKL)
+     // Need a default/reference implementation
+     for (int p = 0; p < batchCount; ++p) {
+       for (int mm = 0; mm < m; ++mm) {
+	 for (int nn = 0; nn < n; ++nn) {
+	   ComplexD c_mn(0.0);
+	   for (int kk = 0; kk < k; ++kk)
+	     c_mn += Amk[mm + kk*lda + p*sda] * Bkn[kk + nn*ldb + p*sdb];
+	   Cmn[mm + nn*ldc + p*sdc] =  (alpha)*c_mn + (beta)*Cmn[mm + nn*ldc + p*sdc];
+	 }
+       }
+     }
+#endif
+  }
 
-  template<class CComplex>
   double benchmark(int M, int N, int K, int BATCH)
   {
     int32_t N_A = M*K*BATCH;
     int32_t N_B = K*N*BATCH;
     int32_t N_C = M*N*BATCH;
-    deviceVector<CComplex> A(N_A); acceleratorMemSet(&A[0],0,N_A*sizeof(CComplex));
+    deviceVector<ComplexD> A(N_A); acceleratorMemSet(&A[0],0,N_A*sizeof(ComplexD));
-    deviceVector<CComplex> B(N_B); acceleratorMemSet(&B[0],0,N_B*sizeof(CComplex));
+    deviceVector<ComplexD> B(N_B); acceleratorMemSet(&B[0],0,N_B*sizeof(ComplexD));
-    deviceVector<CComplex> C(N_C); acceleratorMemSet(&C[0],0,N_C*sizeof(CComplex));
+    deviceVector<ComplexD> C(N_C); acceleratorMemSet(&C[0],0,N_C*sizeof(ComplexD));
-    CComplex alpha(1.0);
+    ComplexD alpha(1.0);
-    CComplex beta (1.0);
+    ComplexD beta (1.0);
     RealD flops = 8.0*M*N*K*BATCH;
-    int ncall=1000;
+    int ncall=10;
-    deviceVector<CComplex *> As(BATCH);
-    deviceVector<CComplex *> Bs(BATCH);
-    deviceVector<CComplex *> Cs(BATCH);
-    for(int b = 0 ; b < BATCH;b++) {
-      CComplex *ptr;
-      ptr = &A[b*M*K];      acceleratorPut(As[b],ptr);
-      ptr = &B[b*K*N];      acceleratorPut(Bs[b],ptr);
-      ptr = &C[b*M*N];      acceleratorPut(Cs[b],ptr);
-    }
-
-    // Warm up call
-    gemmBatched(M,N,K,
-		alpha,
-		As, // m x k 
-		Bs, // k x n
-		beta, 
-		Cs);
-    synchronise();
-
     RealD t0 = usecond();
     for(int i=0;i<ncall;i++){
-      gemmBatched(M,N,K,
+      gemmStridedBatched(M,N,K,
-		  alpha,
+			 alpha,
-		  As, // m x k 
+			 &A[0], // m x k 
-		  Bs, // k x n
+			 &B[0], // k x n
-		  beta, 
+			 beta, 
-		  Cs);
+			 &C[0], // m x n
-      synchronise();
+			 BATCH);
     }
+    synchronise();
     RealD t1 = usecond();
-    RealD bytes = 1.0*sizeof(CComplex)*(M*N*2+N*K+M*K)*BATCH;
+    RealD bytes = 1.0*sizeof(ComplexD)*(M*N*2+N*K+M*K)*BATCH;
     flops = 8.0*M*N*K*BATCH*ncall;
     flops = flops/(t1-t0)/1.e3;
     return flops; // Returns gigaflops
   }
 
+
+
+
 };
 
 NAMESPACE_END(Grid);

Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h

@@ -279,11 +279,11 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
       diagonalize(eval2,lmd2,lme2,Nu,Nm,Nm,Qt,grid);
       _sort.push(eval2,Nm);
-      Glog << "#Ritz value before shift: "<< std::endl;
+      //      Glog << "#Ritz value before shift: "<< std::endl;
       for(int i=0; i<Nm; ++i){
-	std::cout.precision(13);
+	//        std::cout.precision(13);
-	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	//        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 << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
       }
       
       //----------------------------------------------------------------------
@@ -297,8 +297,7 @@ public:
         
         unpackHermitBlockTriDiagMatToEigen(lmd,lme,Nu,Nblock_m,Nm,Nm,BTDM);
 
-        for(int ip=Nk; ip<Nm; ++ip){
+        for(int ip=Nk; ip<Nm; ++ip){ 
-	  Glog << " ip "<<ip<<" / "<<Nm<<std::endl;
           shiftedQRDecompEigen(BTDM,Nu,Nm,eval2[ip],Q);
         }
         
@@ -326,7 +325,7 @@ public:
         Qt = Eigen::MatrixXcd::Identity(Nm,Nm);
         diagonalize(eval2,lmd2,lme2,Nu,Nk,Nm,Qt,grid);
         _sort.push(eval2,Nk);
-	Glog << "#Ritz value after shift: "<< std::endl;
+	//        Glog << "#Ritz value after shift: "<< std::endl;
         for(int i=0; i<Nk; ++i){
 	  //          std::cout.precision(13);
 	  //          std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
@@ -468,10 +467,10 @@ public:
   
     // set initial vector
     for (int i=0; i<Nu; ++i) {
-      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
+      //      Glog << "norm2(src[" << i << "])= "<< norm2(src[i]) << std::endl;
       evec[i] = src[i];
       orthogonalize(evec[i],evec,i);
-      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
+      //      Glog << "norm2(evec[" << i << "])= "<< norm2(evec[i]) << std::endl;
     }
 //    exit(-43);
     
@@ -507,11 +506,11 @@ public:
       Qt = Eigen::MatrixXcd::Identity(Nr,Nr);
       diagonalize(eval2,lmd2,lme2,Nu,Nr,Nr,Qt,grid);
       _sort.push(eval2,Nr);
-      Glog << "#Ritz value: "<< std::endl;
+      //      Glog << "#Ritz value: "<< std::endl;
       for(int i=0; i<Nr; ++i){
-	std::cout.precision(13);
+	//        std::cout.precision(13);
-	std::cout << "[" << std::setw(4)<< std::setiosflags(std::ios_base::right) <<i<<"] ";
+	//        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 << "Rval = "<<std::setw(20)<< std::setiosflags(std::ios_base::left)<< eval2[i] << std::endl;
       }
       
       // Convergence test
@@ -571,7 +570,6 @@ public:
       Glog << fname + " NOT converged ; Summary :\n";
     } else {
       Glog << fname + " CONVERGED ; Summary :\n";
-      Nstop = Nconv_guess; // Just take them all
       // Sort convered eigenpairs.
       std::vector<Field>  Btmp(Nstop,grid); // waste of space replicating
 
@@ -644,7 +642,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 +676,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 +689,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 +714,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) {
@@ -781,7 +779,7 @@ private:
     
     for ( int u=0; u<Nu; ++u ) {
       for (int k=0; k<Nk; ++k ) {
-	//	Glog << "lmd "<<u<<" "<<k<<" "<<lmd[u][k] -conjugate(lmd[u][k])<<std::endl;
+//        Glog << "lmd "<<u<<" "<<k<<" "<<lmd[u][k] -conjugate(lmd[u][k])<<std::endl;
         BlockTriDiag(k,u+(k/Nu)*Nu) = lmd[u][k];
       }
     }
@@ -935,7 +933,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 +951,7 @@ if (1){
         M(u+(k/Nu)*Nu,k-Nu) = lme[u][k-Nu];
       }
     }
-    Glog << "unpackHermitBlockTriDiagMatToEigen() end" << std::endl; 
+    //Glog << "unpackHermitBlockTriDiagMatToEigen() end" << endl; 
   }
  
 
@@ -963,7 +961,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 +977,7 @@ if (1){
         lme[u][k-Nu] = M(u+(k/Nu)*Nu,k-Nu);
       }
     }
-    Glog << "packHermitBlockTriDiagMatfromEigen() end" <<std::endl; 
+    //Glog << "packHermitBlockTriDiagMatfromEigen() end" << endl; 
   }
 
 
@@ -988,7 +986,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 +1002,6 @@ if (1){
                         // lower triangular part used to represent series
                         // of Q sequence.
 
-    Glog << "shiftedQRDecompEigen() Housholder & QR" << '\n'; 
     // equivalent operation of Qprod *= Q
     //M = Eigen::MatrixXcd::Zero(Nm,Nm);
     
@@ -1025,7 +1022,6 @@ if (1){
     
     Mtmp = Eigen::MatrixXcd::Zero(Nm,Nm);
 
-    Glog << "shiftedQRDecompEigen() Mtmp create" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=0; j<Nm-(Nu+1); ++j) {
         for (int k=0; k<Nu+1+j; ++k) {
@@ -1033,7 +1029,6 @@ if (1){
         }
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop1" << '\n'; 
     for (int i=0; i<Nm; ++i) {
       for (int j=Nm-(Nu+1); j<Nm; ++j) {
         for (int k=0; k<Nm; ++k) {
@@ -1041,7 +1036,6 @@ if (1){
         }
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop2" << '\n'; 
     
     //static int ntimes = 2;
     //for (int j=0; j<Nm-(ntimes*Nu); ++j) {
@@ -1067,13 +1061,11 @@ if (1){
         Mtmp(j,i) = conj(Mtmp(i,j));
       }
     }
-    Glog << "shiftedQRDecompEigen() Mtmp loop3" << '\n'; 
 
     for (int i=0; i<Nm; ++i) {
       Mtmp(i,i) = real(Mtmp(i,i)) + Dsh;
     }
     
-    Glog << "shiftedQRDecompEigen() Mtmp loop4" << '\n'; 
     M = Mtmp;
 
     //M = Q.adjoint()*(M*Q);
@@ -1085,7 +1077,7 @@ if (1){
     //  }
     //}
     
-    Glog << "shiftedQRDecompEigen() end" <<std::endl; 
+    //Glog << "shiftedQRDecompEigen() end" << endl; 
   }
 
   void exampleQRDecompEigen(void)

Grid/allocator/MemoryManager.cc

@@ -35,7 +35,7 @@ uint64_t total_host;;
 void MemoryManager::DisplayMallinfo(void)
 {
 #ifdef __linux__
-  struct mallinfo mi; // really want mallinfo2, but glibc version isn't uniform
+  struct mallinfo mi;
   
   mi = mallinfo();
 

Grid/cartesian/Cartesian_base.h

@@ -91,7 +91,6 @@ public:
   ////////////////////////////////////////////////////////////////
   virtual int CheckerBoarded(int dim)=0;
   virtual int CheckerBoard(const Coordinate &site)=0;
-  virtual int CheckerDim(void){ return 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;
   virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0;

Grid/cartesian/Cartesian_red_black.h

@@ -60,7 +60,6 @@ public:
   int              _checker_dim;
   std::vector<int> _checker_board;
 
-  virtual int CheckerDim(void){ return _checker_dim; };
   virtual int CheckerBoarded(int dim){
     if( dim==_checker_dim) return 1;
     else return 0;

Grid/lattice/Lattice_reduction.h

@@ -264,8 +264,24 @@ inline ComplexD rankInnerProduct(const Lattice<vobj> &left,const Lattice<vobj> &
   const uint64_t sites = grid->oSites();
   
   // Might make all code paths go this way.
+#if 0
+  typedef decltype(innerProductD(vobj(),vobj())) inner_t;
+  Vector<inner_t> inner_tmp(sites);
+  auto inner_tmp_v = &inner_tmp[0];
+  {
+    autoView( left_v , left, AcceleratorRead);
+    autoView( right_v,right, AcceleratorRead);
+    // This code could read coalesce
+    // GPU - SIMT lane compliance...
+    accelerator_for( ss, sites, nsimd,{
+	auto x_l = left_v(ss);
+	auto y_l = right_v(ss);
+	coalescedWrite(inner_tmp_v[ss],innerProductD(x_l,y_l));
+    });
+  }
+#else
   typedef decltype(innerProduct(vobj(),vobj())) inner_t;
-  deviceVector<inner_t> inner_tmp(sites);
+  Vector<inner_t> inner_tmp(sites);
   auto inner_tmp_v = &inner_tmp[0];
     
   {
@@ -279,6 +295,7 @@ inline ComplexD rankInnerProduct(const Lattice<vobj> &left,const Lattice<vobj> &
 	coalescedWrite(inner_tmp_v[ss],innerProduct(x_l,y_l));
     });
   }
+#endif
   // This is in single precision and fails some tests
   auto anrm = sumD(inner_tmp_v,sites);  
   nrm = anrm;

Grid/lattice/Lattice_transfer.h

@@ -42,21 +42,50 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
     assert((fine->_rdimensions[d] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[d]); 
   }
 }
+
  
 ////////////////////////////////////////////////////////////////////////////////////////////
 // remove and insert a half checkerboard
 ////////////////////////////////////////////////////////////////////////////////////////////
-
 template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full)
 {
-  acceleratorPickCheckerboard(cb,half,full);
+  half.Checkerboard() = cb;
+
+  autoView( half_v, half, CpuWrite);
+  autoView( full_v, full, CpuRead);
+  thread_for(ss, full.Grid()->oSites(),{
+    int cbos;
+    Coordinate coor;
+    full.Grid()->oCoorFromOindex(coor,ss);
+    cbos=half.Grid()->CheckerBoard(coor);
+
+    if (cbos==cb) {
+      int ssh=half.Grid()->oIndex(coor);
+      half_v[ssh] = full_v[ss];
+    }
+  });
 }
 template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half)
 {
-  acceleratorSetCheckerboard(full,half);
+  int cb = half.Checkerboard();
+  autoView( half_v , half, CpuRead);
+  autoView( full_v , full, CpuWrite);
+  thread_for(ss,full.Grid()->oSites(),{
+
+    Coordinate coor;
+    int cbos;
+
+    full.Grid()->oCoorFromOindex(coor,ss);
+    cbos=half.Grid()->CheckerBoard(coor);
+      
+    if (cbos==cb) {
+      int ssh=half.Grid()->oIndex(coor);
+      full_v[ss]=half_v[ssh];
+    }
+  });
 }
 
-template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int dummy=0)
+template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full, int checker_dim_half=0)
 {
   half.Checkerboard() = cb;
   autoView(half_v, half, AcceleratorWrite);
@@ -66,7 +95,6 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
-  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss, full.Grid()->oSites(),full.Grid()->Nsimd(),{
     
     Coordinate coor;
@@ -91,7 +119,7 @@ template<class vobj> inline void acceleratorPickCheckerboard(int cb,Lattice<vobj
     }
   });
 }
-template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int dummy=0)
+template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half, int checker_dim_half=0)
 {
   int cb = half.Checkerboard();
   autoView(half_v , half, AcceleratorRead);
@@ -101,7 +129,6 @@ template<class vobj> inline void acceleratorSetCheckerboard(Lattice<vobj> &full,
   unsigned long ndim_half          = half.Grid()->_ndimension;
   Coordinate checker_dim_mask_half = half.Grid()->_checker_dim_mask;
   Coordinate ostride_half          = half.Grid()->_ostride;
-  int checker_dim_half             = half.Grid()->CheckerDim();
   accelerator_for(ss,full.Grid()->oSites(),full.Grid()->Nsimd(),{
 
     Coordinate coor;

Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -86,8 +86,13 @@ public:
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);
 
+#if defined(GRID_CUDA) || defined(GRID_HIP)  || defined(GRID_SYCL)
+    acceleratorSetCheckerboard(Force,ForceE);
+    acceleratorSetCheckerboard(Force,ForceO);
+#else
     setCheckerboard(Force,ForceE); 
     setCheckerboard(Force,ForceO);
+#endif
     Force=-Force;
 
     delete forcecb;
@@ -130,8 +135,13 @@ public:
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);
 
+#if defined(GRID_CUDA) || defined(GRID_HIP)  || defined(GRID_SYCL)
+    acceleratorSetCheckerboard(Force,ForceE);
+    acceleratorSetCheckerboard(Force,ForceO);
+#else
     setCheckerboard(Force,ForceE); 
     setCheckerboard(Force,ForceO);
+#endif
     Force=-Force;
 
     delete forcecb;

Grid/qcd/smearing/GaugeConfigurationMasked.h

@@ -32,9 +32,7 @@ private:
   //  Smear_Stout<Gimpl> *StoutSmearing;
   //  std::vector<GaugeField> SmearedSet;
   
-  GridRedBlackCartesian * UrbGrid; // keep a copy of the redblack grid for life of object
   std::vector<LatticeLorentzComplex> masks;
-  std::vector<int> cbs;
 
   typedef typename SU3Adjoint::AMatrix AdjMatrix;
   typedef typename SU3Adjoint::LatticeAdjMatrix  AdjMatrixField;
@@ -149,25 +147,6 @@ private:
     }
     pokeLorentz(Fdet, Fdet_pol, nu);
   }
-
-  void Compute_MpInvJx_dNxxdSy(int cb,
-			       const GaugeLinkField &PlaqL,
-			       const GaugeLinkField &PlaqR,
-			       AdjMatrixField MpInvJx,
-			       AdjVectorField &Fdet2 )
-  {
-    GaugeLinkField PlaqLeo(UrbGrid);
-    GaugeLinkField PlaqReo(UrbGrid);
-    AdjMatrixField MpInvJxeo(UrbGrid);
-    AdjVectorField Fdet2eo(UrbGrid);
-    pickCheckerboard(cb,PlaqLeo,PlaqL);
-    pickCheckerboard(cb,PlaqReo,PlaqR);
-    pickCheckerboard(cb,MpInvJxeo,MpInvJx);
-    Fdet2eo.Checkerboard()=cb;
-    Compute_MpInvJx_dNxxdSy(PlaqLeo,PlaqReo,MpInvJxeo,Fdet2eo);
-    setCheckerboard(Fdet2,Fdet2eo);
-  }
-  
   void Compute_MpInvJx_dNxxdSy(const GaugeLinkField &PlaqL,const GaugeLinkField &PlaqR, AdjMatrixField MpInvJx,AdjVectorField &Fdet2 )
   {
     GaugeLinkField UtaU(PlaqL.Grid());
@@ -299,9 +278,8 @@ public:
     ////////////////////////////////////////////////////////////////////////////////
     // Mask the gauge field
     ////////////////////////////////////////////////////////////////////////////////
-    int cb = cbs[smr];
     auto mask=PeekIndex<LorentzIndex>(masks[smr],mu); // the cb mask
-    
+
     Umsk = U;
     ApplyMask(Umsk,smr);
     Utmp = peekLorentz(Umsk,mu);
@@ -464,7 +442,7 @@ public:
     AdjMatrixField MpInvJx_nu(grid);
     MpInvJx = (-1.0)*MpAdInv * JxAd;// rho is on the plaq factor
 
-    Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx,FdetV);
+    Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
     Fdet2_mu=FdetV;
     Fdet1_mu=Zero();
     
@@ -521,7 +499,7 @@ public:
 
 	time=-usecond();
 	PlaqR=(-1.0)*PlaqR;
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx,FdetV);
 	Fdet2_nu = FdetV;
 	time+=usecond();
 	std::cout << GridLogMessage << "Compute_MpInvJx_dNxxSy (occurs 6x) took "<<time<< " us"<<std::endl;
@@ -542,7 +520,7 @@ public:
 	
 
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 	
 	///////////////// -ve nu /////////////////
@@ -561,7 +539,7 @@ public:
 	Fdet1_nu = Fdet1_nu + transpose(Nxy)*dJdXe_nMpInv_y;
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 	
 	// x==
@@ -582,7 +560,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,mu,-1);
 	MpInvJx_nu = Cshift(MpInvJx_nu,nu,1);
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_nu = Fdet2_nu+FdetV;
 
 	/////////////////////////////////////////////////////////////////////
@@ -611,7 +589,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,-1);
 
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
 
 	//  __
@@ -631,7 +609,7 @@ public:
 
 	MpInvJx_nu = Cshift(MpInvJx,nu,1);
 
-	Compute_MpInvJx_dNxxdSy(cb,PlaqL,PlaqR,MpInvJx_nu,FdetV);
+	Compute_MpInvJx_dNxxdSy(PlaqL,PlaqR,MpInvJx_nu,FdetV);
 	Fdet2_mu = Fdet2_mu+FdetV;
 	
       }
@@ -953,10 +931,6 @@ private:
 public:
 
   /* Standard constructor */
-  virtual ~SmearedConfigurationMasked()
-  {
-    delete UrbGrid;
-  }
   SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout)
     : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
   {
@@ -965,6 +939,7 @@ public:
     // was resized in base class
     assert(this->SmearedSet.size()==Nsmear);
     
+    GridRedBlackCartesian * UrbGrid;
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);
     LatticeComplex one(_UGrid); one = ComplexD(1.0,0.0);
     LatticeComplex tmp(_UGrid);
@@ -972,11 +947,10 @@ public:
     for (unsigned int i = 0; i < this->smearingLevels; ++i) {
 
       masks.push_back(*(new LatticeLorentzComplex(_UGrid)));
+
       int mu= (i/2) %Nd;
       int cb= (i%2);
       LatticeComplex tmpcb(UrbGrid);
-
-      cbs.push_back(cb);
 	
       masks[i]=Zero();
       ////////////////////
@@ -988,6 +962,7 @@ public:
       PokeIndex<LorentzIndex>(masks[i],tmp, mu);
 	
     }
+    delete UrbGrid;
   }
   
   virtual void smeared_force(GaugeField &SigmaTilde) 

Grid/qcd/utils/GaugeGroup.h

@@ -418,32 +418,32 @@ static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in,
   int hNNm1= NNm1/2;
   RealD sqrt_2 = sqrt(2.0);
   Complex ci(0.0,1.0);
-
+  for(int su2Index=0;su2Index<hNNm1;su2Index++){
-  const int nsimd=  Matrix::Nsimd();
+    int i1, i2;
-  accelerator_for(ss,grid->oSites(),nsimd,{
+    su2SubGroupIndex(i1, i2, su2Index);
-      for(int su2Index=0;su2Index<hNNm1;su2Index++){
+    int ax = su2Index*2;
-	int i1, i2;
+    int ay = su2Index*2+1;
-	su2SubGroupIndex(i1, i2, su2Index);
+    accelerator_for(ss,grid->oSites(),1,{
-	int ax = su2Index*2;
-	int ay = su2Index*2+1;
 	// in is traceless ANTI-hermitian whereas Grid generators are Hermitian.
 	// trace( Ta x Ci in)
 	// Bet I need to move to real part with mult by -i
-	coalescedWrite(out_v[ss]()()(ax,b),0.5*(real(in_v(ss)()()(i2,i1)) - real(in_v(ss)()()(i1,i2))));
+	out_v[ss]()()(ax,b) = 0.5*(real(in_v[ss]()()(i2,i1)) - real(in_v[ss]()()(i1,i2)));
-	coalescedWrite(out_v[ss]()()(ay,b),0.5*(imag(in_v(ss)()()(i1,i2)) + imag(in_v(ss)()()(i2,i1))));
+	out_v[ss]()()(ay,b) = 0.5*(imag(in_v[ss]()()(i1,i2)) + imag(in_v[ss]()()(i2,i1)));
-      }
+      });
-      for(int diagIndex=0;diagIndex<N-1;diagIndex++){
+  }
-	int k = diagIndex + 1; // diagIndex starts from 0
+  for(int diagIndex=0;diagIndex<N-1;diagIndex++){
-	int a = NNm1+diagIndex;
+    int k = diagIndex + 1; // diagIndex starts from 0
-	RealD scale = 1.0/sqrt(2.0*k*(k+1));
+    int a = NNm1+diagIndex;
-	auto tmp = in_v(ss)()()(0,0);
+    RealD scale = 1.0/sqrt(2.0*k*(k+1));
+    accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
+	auto tmp = in_v[ss]()()(0,0);
 	for(int i=1;i<k;i++){
-	  tmp=tmp+in_v(ss)()()(i,i);
+	  tmp=tmp+in_v[ss]()()(i,i);
 	}
-	tmp = tmp - in_v(ss)()()(k,k)*k;
+	tmp = tmp - in_v[ss]()()(k,k)*k;
-	coalescedWrite(out_v[ss]()()(a,b),imag(tmp) * scale);
+	out_v[ss]()()(a,b) =imag(tmp) * scale;
-      }
+      });
-    });
+    }
 }
 
   

Grid/qcd/utils/SUn.impl.h

@@ -118,7 +118,7 @@ static void generatorDiagonal(int diagIndex, iGroupMatrix<cplx> &ta) {
 ////////////////////////////////////////////////////////////////////////
 // Map a su2 subgroup number to the pair of rows that are non zero
 ////////////////////////////////////////////////////////////////////////
-static accelerator_inline void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
+static void su2SubGroupIndex(int &i1, int &i2, int su2_index, GroupName::SU) {
   assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2));
 
   int spare = su2_index;

Grid/tensors/Tensor_class.h

@@ -460,9 +460,3 @@ void vprefetch(const iMatrix<v, N> &vv) {
 
 NAMESPACE_END(Grid);
 
-
-#ifdef GRID_SYCL
-template<class vec> struct sycl::is_device_copyable<Grid::iScalar<vec> > : public std::true_type {};
-template<class vec,int N> struct sycl::is_device_copyable<Grid::iVector<vec,N> > : public std::true_type {};
-template<class vec,int N> struct sycl::is_device_copyable<Grid::iMatrix<vec,N> > : public std::true_type {};
-#endif

benchmarks/Benchmark_usqcd.cc

@@ -261,25 +261,23 @@ public:
     fprintf(FP,"\n\n");
   };
 
-  template<class CComplex>
+
   static void BLAS(void)
   {
     //int nbasis, int nrhs, int coarseVol
     int  basis[] = { 16,32,64 };
-    int  rhs[]   = { 8,12,16 };
+    int  rhs[]   = { 8,16,32 };
-    int  vol  = 8*8*8*8;
+    int  vol  = 4*4*4*4;
-    int  blk  = 4*4*4*4;
 
     GridBLAS blas;
-
+    
-    int fpbits = sizeof(CComplex)*4;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "= batched GEMM fp"<<fpbits<<std::endl;
+    std::cout<<GridLogMessage << "= batched GEMM (double precision) "<<std::endl;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     std::cout<<GridLogMessage << "  M  "<<"\t\t"<<"N"<<"\t\t\t"<<"K"<<"\t\t"<<"Gflop/s / rank (coarse mrhs)"<<std::endl;
     std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
   
-    fprintf(FP,"GEMM\n\n M, N, K, BATCH, GF/s per rank fp%d\n",fpbits);
+    fprintf(FP,"GEMM\n\n M, N, K, BATCH, GF/s per rank\n");
 
     for(int b=0;b<3;b++){
     for(int r=0;r<3;r++){
@@ -287,7 +285,7 @@ public:
       int N=rhs[r];
       int K=basis[b];
       int BATCH=vol;
-      double p=blas.benchmark<CComplex>(M,N,K,BATCH);
+      double p=blas.benchmark(M,N,K,BATCH);
 
       fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p);
       
@@ -301,9 +299,9 @@ public:
     for(int r=0;r<3;r++){
       int M=basis[b];
       int N=rhs[r];
-      int K=blk;
+      int K=vol;
       int BATCH=vol;
-      double p=blas.benchmark<CComplex>(M,N,K,BATCH);
+      double p=blas.benchmark(M,N,K,BATCH);
 
       fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p);
       std::cout<<GridLogMessage<<std::setprecision(3) 
@@ -315,10 +313,10 @@ public:
     for(int b=0;b<3;b++){
     for(int r=0;r<3;r++){
       int M=rhs[r];
-      int N=blk;
+      int N=vol;
       int K=basis[b];
       int BATCH=vol;
-      double p=blas.benchmark<CComplex>(M,N,K,BATCH);
+      double p=blas.benchmark(M,N,K,BATCH);
 
       fprintf(FP,"%d, %d, %d, %d, %f\n", M, N, K, BATCH, p);
       std::cout<<GridLogMessage<<std::setprecision(3) 
@@ -869,7 +867,6 @@ int main (int argc, char ** argv)
   int do_memory=1;
   int do_comms =1;
   int do_blas  =1;
-  int do_dslash=1;
 
   int sel=4;
   std::vector<int> L_list({8,12,16,24,32});
@@ -880,7 +877,6 @@ int main (int argc, char ** argv)
   std::vector<double> staggered;
 
   int Ls=1;
-  if (do_dslash){
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
   std::cout<<GridLogMessage << " Clover dslash 4D vectorised (temporarily Wilson)" <<std::endl;
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
@@ -905,7 +901,6 @@ int main (int argc, char ** argv)
     staggered.push_back(result);
   }
 
-
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
   std::cout<<GridLogMessage << " Summary table Ls="<<Ls <<std::endl;
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
@@ -914,33 +909,8 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]<<" \t\t "<<dwf4[l] << " \t\t "<< staggered[l]<<std::endl;
   }
   std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  }
 
   int NN=NN_global;
-  if(do_dslash){
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Per Node Summary table Ls="<<Ls <<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " L \t\t Clover\t\t DWF4\t\t Staggered (GF/s per node)" <<std::endl;
-    fprintf(FP,"Per node summary table\n");
-    fprintf(FP,"\n");
-    fprintf(FP,"L , Wilson, DWF4, Staggered, GF/s per node\n");
-    fprintf(FP,"\n");
-    for(int l=0;l<L_list.size();l++){
-      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]/NN<<" \t "<<dwf4[l]/NN<< " \t "<<staggered[l]/NN<<std::endl;
-      fprintf(FP,"%d , %.0f, %.0f, %.0f\n",L_list[l],clover[l]/NN/1000.,dwf4[l]/NN/1000.,staggered[l]/NN/1000.);
-    }
-    fprintf(FP,"\n");
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Comparison point     result: "  << 0.5*(dwf4[sel]+dwf4[selm1])/NN << " Mflop/s per node"<<std::endl;
-    std::cout<<GridLogMessage << " Comparison point is 0.5*("<<dwf4[sel]/NN<<"+"<<dwf4[selm1]/NN << ") "<<std::endl;
-    std::cout<<std::setprecision(3);
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  }
-
-  
   if ( do_memory ) {
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     std::cout<<GridLogMessage << " Memory benchmark " <<std::endl;
@@ -948,6 +918,15 @@ int main (int argc, char ** argv)
     Benchmark::Memory();
   }
 
+  if ( do_blas ) {
+#if defined(GRID_CUDA) || defined(GRID_HIP)     || defined(GRID_SYCL)   
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+    std::cout<<GridLogMessage << " Batched BLAS benchmark " <<std::endl;
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+    Benchmark::BLAS();
+#endif
+  }
+
   if ( do_su4 ) {
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
     std::cout<<GridLogMessage << " SU(4) benchmark " <<std::endl;
@@ -962,14 +941,28 @@ int main (int argc, char ** argv)
     Benchmark::Comms();
   }
 
-  if ( do_blas ) {
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Batched BLAS benchmark " <<std::endl;
+    std::cout<<GridLogMessage << " Per Node Summary table Ls="<<Ls <<std::endl;
     std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    Benchmark::BLAS<ComplexD>();
+    std::cout<<GridLogMessage << " L \t\t Clover\t\t DWF4\t\t Staggered (GF/s per node)" <<std::endl;
-    Benchmark::BLAS<ComplexF>();
+    fprintf(FP,"Per node summary table\n");
-  }
+    fprintf(FP,"\n");
-  
+    fprintf(FP,"L , Wilson, DWF4, Staggered, GF/s per node\n");
+    fprintf(FP,"\n");
+    for(int l=0;l<L_list.size();l++){
+      std::cout<<GridLogMessage << L_list[l] <<" \t\t "<< clover[l]/NN<<" \t "<<dwf4[l]/NN<< " \t "<<staggered[l]/NN<<std::endl;
+      fprintf(FP,"%d , %.0f, %.0f, %.0f\n",L_list[l],clover[l]/NN/1000.,dwf4[l]/NN/1000.,staggered[l]/NN/1000.);
+    }
+    fprintf(FP,"\n");
+
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+    std::cout<<GridLogMessage << " Comparison point     result: "  << 0.5*(dwf4[sel]+dwf4[selm1])/NN << " Mflop/s per node"<<std::endl;
+    std::cout<<GridLogMessage << " Comparison point is 0.5*("<<dwf4[sel]/NN<<"+"<<dwf4[selm1]/NN << ") "<<std::endl;
+    std::cout<<std::setprecision(3);
+    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
+
   Grid_finalize();
   fclose(FP);
 }

systems/Aurora/config-command-leak

@@ -1,23 +0,0 @@
-source ~/spack/share/spack/setup-env.sh 
-spack load c-lime
-export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
-export TCMALLOC=`spack find --paths gperftools | grep ^gperftools | awk '{print $2}' `
-export LD_LIBRARY_PATH=${TCMALLOC}/lib:$LD_LIBRARY_PATH
-
-../../configure \
-	--enable-debug \
-	--enable-simd=GPU \
-	--enable-gen-simd-width=64 \
-	--enable-comms=mpi-auto \
-	--disable-gparity \
-	--disable-fermion-reps \
-	--with-lime=$CLIME \
-	--enable-shm=nvlink \
-	--enable-accelerator=sycl \
-	--enable-accelerator-aware-mpi=yes\
-	--enable-unified=no \
-	MPICXX=mpicxx \
-	CXX=icpx \
-	LDFLAGS="-fiopenmp -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel -fsycl  -lsycl -Xarch_host -fsanitize=leak -fsycl-device-code-split=per_kernel" \
-	CXXFLAGS="-fiopenmp -fsycl-unnamed-lambda -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -qmkl=parallel -Xarch_host  -fsycl -fsanitize=leak "
-

systems/Aurora/sourceme.sh

@@ -1,25 +1,13 @@
 source ~/spack/share/spack/setup-env.sh 
 spack load c-lime
+
 export CLIME=`spack find --paths c-lime | grep ^c-lime | awk '{print $2}' `
-#spack load libefence
+#export LD_LIBRARY_PATH=${TCMALLOC}/lib:$LD_LIBRARY_PATH
-#export EFENCE=`spack find --paths libefence | grep ^libefence | awk '{print $2}' `
+
-#export LD_LIBRARY_PATH=${EFENCE}/lib:$LD_LIBRARY_PATH
-#spack load gperftools
-export TCMALLOC=/home/paboyle/gperftools/install
-export LD_LIBRARY_PATH=${TCMALLOC}/lib:$LD_LIBRARY_PATH
 export INTELGT_AUTO_ATTACH_DISABLE=1
 
 #export ONEAPI_DEVICE_SELECTOR=level_zero:0.0
-#module load oneapi/release/2023.12.15.001
-#module use /soft/modulefiles
-#module load intel_compute_runtime/release/agama-devel-682.22
 
-#export FI_CXI_DEFAULT_CQ_SIZE=131072
-#export FI_CXI_CQ_FILL_PERCENT=20
-#export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
-#export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-intel-enable-auto-large-GRF-mode"
-
-#
 # -ftarget-register-alloc-mode=pvc:default 
 # -ftarget-register-alloc-mode=pvc:small
 # -ftarget-register-alloc-mode=pvc:large
@@ -32,9 +20,4 @@ export http_proxy=http://proxy.alcf.anl.gov:3128
 export https_proxy=http://proxy.alcf.anl.gov:3128
 git config --global http.proxy http://proxy.alcf.anl.gov:3128
 
-#source ~/spack/share/spack/setup-env.sh
-#spack load gperftools
-#export TCMALLOC=`spack find --paths gperftools | grep ^gperftools | awk '{print $2}' `
-#export LD_LIBRARY_PATH=${TCMALLOC}/lib:$LD_LIBRARY_PATH
-
 export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"

systems/Frontier/sourceme.sh

@@ -3,7 +3,7 @@ spack load c-lime
 module load emacs 
 module load PrgEnv-gnu
 module load rocm
-module load cray-mpich
+module load cray-mpich/8.1.23
 module load gmp
 module load cray-fftw
 module load craype-accel-amd-gfx90a

tests/debug/Test_8888.cc

@@ -1,118 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_general_coarse_hdcg.cc
-
-    Copyright (C) 2023
-
-Author: Peter Boyle <pboyle@bnl.gov>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#include <Grid/Grid.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczos.h>
-#include <Grid/algorithms/iterative/ImplicitlyRestartedBlockLanczosCoarse.h>
-#include <Grid/algorithms/iterative/AdefMrhs.h>
-
-using namespace std;
-using namespace Grid;
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  const int Ls=8;
-  const int nbasis = 40;
-  const int cb = 0 ;
-  RealD mass=0.01;
-  RealD M5=1.8;
-  RealD b=1.0;
-  RealD c=0.0;
-
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
-								   GridDefaultSimd(Nd,vComplex::Nsimd()),
-								   GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  ///////////////////////// RNGs /////////////////////////////////
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> seeds5({5,6,7,8});
-  std::vector<int> cseeds({5,6,7,8});
-
-  GridParallelRNG          RNG5(FGrid);   RNG5.SeedFixedIntegers(seeds5);
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-
-  ///////////////////////// Configuration /////////////////////////////////
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_EODWF_lat.125");
-  NerscIO::readConfiguration(Umu,header,file);
-
-  //////////////////////// Fermion action //////////////////////////////////
-  MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
-
-  MdagMLinearOperator<MobiusFermionD, LatticeFermion> HermOp(Ddwf);
-
-  
-  std::cout << "**************************************"<<std::endl;
-  std::cout << "         Fine Power method            "<<std::endl;
-  std::cout << "**************************************"<<std::endl;
-
-  LatticeFermionD pm_src(FGrid);
-  pm_src = ComplexD(1.0);
-  PowerMethod<LatticeFermionD>       fPM;
-  fPM(HermOp,pm_src);
-
-  
-  std::cout << "**************************************"<<std::endl;
-  std::cout << "         Fine Lanczos  (poly, low)    "<<std::endl;
-  std::cout << "**************************************"<<std::endl;
-  
-  int Nk=80;
-  int Nm=Nk*3;
-  int Nstop=8;
-  int Nconv_test_interval=1;
-  
-  //  Chebyshev<LatticeFermionD>      IRLChebyLo(0.2,64.0,201);  // 1 iter
-  Chebyshev<LatticeFermionD>      IRLChebyLo(0.0,55.0,101);  // 1 iter
-  FunctionHermOp<LatticeFermionD>    PolyOp(IRLChebyLo,HermOp);
-  PlainHermOp<LatticeFermionD>          Op(HermOp);
-
-  ImplicitlyRestartedLanczos IRL(PolyOp,
-				 Op,
-				 Nk, // sought vecs
-				 Nk, // sought vecs
-				 Nm, // spare vecs
-				 1.0e-8,
-				 10 // Max iterations
-				 );
-
-  int Nconv;
-  std::vector<RealD>            eval(Nm);
-  std::vector<LatticeFermionD>     evec(Nm,FGrid);
-  LatticeFermionD     irl_src(FGrid);
-
-  IRL.calc(eval,evec,irl_src,Nconv);
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse.cc

@@ -244,7 +244,7 @@ int main (int argc, char ** argv)
 
   GridCartesian *CoarseMrhs = new GridCartesian(rhLatt,rhSimd,rhMpi); 
 
-#if 0  
+  
   MultiGeneralCoarsenedMatrix mrhs(LittleDiracOp,CoarseMrhs);
   typedef decltype(mrhs) MultiGeneralCoarsenedMatrix_t;
   
@@ -307,8 +307,7 @@ int main (int argc, char ** argv)
     rh_res= Zero();
     mrhsCG(MrhsCoarseOp,rh_src,rh_res);
   }
-
+  
-#endif
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<std::endl;
   std::cout<<GridLogMessage<<"*******************************************"<<std::endl;

tests/debug/Test_general_coarse_hdcg_phys48_mixed.cc

@@ -145,7 +145,7 @@ int main (int argc, char ** argv)
   Grid_init(&argc,&argv);
 
   const int Ls=24;
-  const int nbasis = 62;
+  const int nbasis = 60;
   const int cb = 0 ;
   RealD mass=0.00078;
   RealD M5=1.8;
@@ -160,7 +160,7 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
   // Construct a coarsened grid with 4^4 cell
-  Coordinate Block({4,4,6,4});
+  Coordinate Block({4,4,4,4});
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/Block[d];

tests/debug/Test_general_coarse_hdcg_phys96_mixed.cc

@@ -160,8 +160,7 @@ int main (int argc, char ** argv)
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
   // Construct a coarsened grid with 4^4 cell
-  //  Coordinate Block({4,4,6,4});
+  Coordinate Block({4,4,6,6});
-  Coordinate Block({4,4,4,4});
   Coordinate clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/Block[d];
@@ -218,7 +217,7 @@ int main (int argc, char ** argv)
   std::string evec_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/evecs.scidac");
   std::string eval_file("/lustre/orion/phy157/proj-shared/phy157_dwf/paboyle/eval.xml");
   bool load_agg=true;
-  bool load_refine=true;
+  bool load_refine=false;
   bool load_mat=false;
   bool load_evec=false;
 
@@ -277,25 +276,17 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
 
   typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
-  //  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,301);  // 4.4.6.4
+  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,301);  // 1 iter
-  //  Chebyshev<CoarseVector>      IRLCheby(0.0012,42.0,501);  // for 4.4.4.4 blocking 350 evs
-  //  Chebyshev<CoarseVector>      IRLCheby(0.0014,42.0,501);  // for 4.4.4.4 blocking 700 evs
-  //  Chebyshev<CoarseVector>      IRLCheby(0.002,42.0,501);  // for 4.4.4.4 blocking 1226 evs
-  //  Chebyshev<CoarseVector>      IRLCheby(0.0025,42.0,501);  // for 4.4.4.4 blocking 1059 evs
-							  //							  3e-4,2);
-  Chebyshev<CoarseVector>      IRLCheby(0.0018,42.0,301);  // for 4.4.4.4 blocking  // 790 evs
-  
   MrhsHermMatrix MrhsCoarseOp     (mrhs);
 
   CoarseVector pm_src(CoarseMrhs);
   pm_src = ComplexD(1.0);
   PowerMethod<CoarseVector>       cPM;   cPM(MrhsCoarseOp,pm_src);
 
-  //  int Nk=nrhs*30; // 4.4.6.4
+  int Nk=nrhs*30;
   //  int Nk=nrhs*80;
-  int Nk=nrhs*60; // 720
+  int Nm=Nk*4;
-  int Nm=Nk*4;    // 2880 ; generally finishes at 1440
+  int Nstop=Nk;
-  int Nstop=512;
   int Nconv_test_interval=1;
   
   ImplicitlyRestartedBlockLanczosCoarse<CoarseVector> IRL(MrhsCoarseOp,
@@ -308,7 +299,7 @@ int main (int argc, char ** argv)
 							  nrhs,
 							  Nk,
 							  Nm,
-							  3e-4,2);
+							  1e-4,20);
 
   std::vector<RealD>            eval(Nm);
   std::vector<CoarseVector>     evec(Nm,Coarse5d);
@@ -340,7 +331,7 @@ int main (int argc, char ** argv)
   // Extra HDCG parameters
   //////////////////////////
   int maxit=3000;
-  ConjugateGradient<CoarseVector>  CG(7.5e-2,maxit,false);
+  ConjugateGradient<CoarseVector>  CG(5.0e-2,maxit,false);
   RealD lo=2.0;
   int ord = 7;