Merge b02d022993 into ee1b8bbdbd

Grid/Grid_Eigen_Dense.h

@@ -34,7 +34,7 @@
 #pragma push_macro("__SYCL_DEVICE_ONLY__")
 #undef __SYCL_DEVICE_ONLY__
 #define EIGEN_DONT_VECTORIZE
-#undef EIGEN_USE_SYCL
+//#undef EIGEN_USE_SYCL
 #define __SYCL__REDEFINE__
 #endif
 

Grid/algorithms/FFT.h

@@ -29,7 +29,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #define _GRID_FFT_H_
 
 #ifdef HAVE_FFTW
-#if defined(USE_MKL) || defined(GRID_SYCL)
+#ifdef USE_MKL
 #include <fftw/fftw3.h>
 #else
 #include <fftw3.h>

Grid/algorithms/approx/Zolotarev.cc

@@ -293,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k,
  * Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and
  * type = 1 for the approximation which is infinite at x = 0. */
 
-zolotarev_data* zolotarev(ZOLO_PRECISION epsilon, int n, int type) {
+zolotarev_data* zolotarev(PRECISION epsilon, int n, int type) {
   INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F,
     l, invlambda, xi, xisq, *tv, s, opl;
   int m, czero, ts;
@@ -375,12 +375,12 @@ zolotarev_data* zolotarev(ZOLO_PRECISION epsilon, int n, int type) {
   construct_partfrac(d);
   construct_contfrac(d);
 
-  /* Converting everything to ZOLO_PRECISION for external use only */
+  /* Converting everything to PRECISION for external use only */
 
   zd = (zolotarev_data*) malloc(sizeof(zolotarev_data));
-  zd -> A = (ZOLO_PRECISION) d -> A;
-  zd -> Delta = (ZOLO_PRECISION) d -> Delta;
-  zd -> epsilon = (ZOLO_PRECISION) d -> epsilon;
+  zd -> A = (PRECISION) d -> A;
+  zd -> Delta = (PRECISION) d -> Delta;
+  zd -> epsilon = (PRECISION) d -> epsilon;
   zd -> n = d -> n;
   zd -> type = d -> type;
   zd -> dn = d -> dn;
@@ -390,24 +390,24 @@ zolotarev_data* zolotarev(ZOLO_PRECISION epsilon, int n, int type) {
   zd -> deg_num = d -> deg_num;
   zd -> deg_denom = d -> deg_denom;
 
-  zd -> a = (ZOLO_PRECISION*) malloc(zd -> dn * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (ZOLO_PRECISION) d -> a[m];
+  zd -> a = (PRECISION*) malloc(zd -> dn * sizeof(PRECISION));
+  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (PRECISION) d -> a[m];
   free(d -> a);
 
-  zd -> ap = (ZOLO_PRECISION*) malloc(zd -> dd * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (ZOLO_PRECISION) d -> ap[m];
+  zd -> ap = (PRECISION*) malloc(zd -> dd * sizeof(PRECISION));
+  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (PRECISION) d -> ap[m];
   free(d -> ap);
 
-  zd -> alpha = (ZOLO_PRECISION*) malloc(zd -> da * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (ZOLO_PRECISION) d -> alpha[m];
+  zd -> alpha = (PRECISION*) malloc(zd -> da * sizeof(PRECISION));
+  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (PRECISION) d -> alpha[m];
   free(d -> alpha);
 
-  zd -> beta = (ZOLO_PRECISION*) malloc(zd -> db * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (ZOLO_PRECISION) d -> beta[m];
+  zd -> beta = (PRECISION*) malloc(zd -> db * sizeof(PRECISION));
+  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (PRECISION) d -> beta[m];
   free(d -> beta);
 
-  zd -> gamma = (ZOLO_PRECISION*) malloc(zd -> n * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (ZOLO_PRECISION) d -> gamma[m];
+  zd -> gamma = (PRECISION*) malloc(zd -> n * sizeof(PRECISION));
+  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (PRECISION) d -> gamma[m];
   free(d -> gamma);
 
   free(d);
@@ -426,7 +426,7 @@ void zolotarev_free(zolotarev_data *zdata)
 }
 
 
-zolotarev_data* higham(ZOLO_PRECISION epsilon, int n) {
+zolotarev_data* higham(PRECISION epsilon, int n) {
   INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq;
   int m, czero;
   zolotarev_data *zd;
@@ -481,9 +481,9 @@ zolotarev_data* higham(ZOLO_PRECISION epsilon, int n) {
   /* Converting everything to PRECISION for external use only */
 
   zd = (zolotarev_data*) malloc(sizeof(zolotarev_data));
-  zd -> A = (ZOLO_PRECISION) d -> A;
-  zd -> Delta = (ZOLO_PRECISION) d -> Delta;
-  zd -> epsilon = (ZOLO_PRECISION) d -> epsilon;
+  zd -> A = (PRECISION) d -> A;
+  zd -> Delta = (PRECISION) d -> Delta;
+  zd -> epsilon = (PRECISION) d -> epsilon;
   zd -> n = d -> n;
   zd -> type = d -> type;
   zd -> dn = d -> dn;
@@ -493,24 +493,24 @@ zolotarev_data* higham(ZOLO_PRECISION epsilon, int n) {
   zd -> deg_num = d -> deg_num;
   zd -> deg_denom = d -> deg_denom;
 
-  zd -> a = (ZOLO_PRECISION*) malloc(zd -> dn * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (ZOLO_PRECISION) d -> a[m];
+  zd -> a = (PRECISION*) malloc(zd -> dn * sizeof(PRECISION));
+  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (PRECISION) d -> a[m];
   free(d -> a);
 
-  zd -> ap = (ZOLO_PRECISION*) malloc(zd -> dd * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (ZOLO_PRECISION) d -> ap[m];
+  zd -> ap = (PRECISION*) malloc(zd -> dd * sizeof(PRECISION));
+  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (PRECISION) d -> ap[m];
   free(d -> ap);
 
-  zd -> alpha = (ZOLO_PRECISION*) malloc(zd -> da * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (ZOLO_PRECISION) d -> alpha[m];
+  zd -> alpha = (PRECISION*) malloc(zd -> da * sizeof(PRECISION));
+  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (PRECISION) d -> alpha[m];
   free(d -> alpha);
 
-  zd -> beta = (ZOLO_PRECISION*) malloc(zd -> db * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (ZOLO_PRECISION) d -> beta[m];
+  zd -> beta = (PRECISION*) malloc(zd -> db * sizeof(PRECISION));
+  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (PRECISION) d -> beta[m];
   free(d -> beta);
 
-  zd -> gamma = (ZOLO_PRECISION*) malloc(zd -> n * sizeof(ZOLO_PRECISION));
-  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (ZOLO_PRECISION) d -> gamma[m];
+  zd -> gamma = (PRECISION*) malloc(zd -> n * sizeof(PRECISION));
+  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (PRECISION) d -> gamma[m];
   free(d -> gamma);
 
   free(d);
@@ -523,17 +523,17 @@ NAMESPACE_END(Grid);
 #ifdef TEST
 
 #undef ZERO
-#define ZERO ((ZOLO_PRECISION) 0)
+#define ZERO ((PRECISION) 0)
 #undef ONE
-#define ONE ((ZOLO_PRECISION) 1)
+#define ONE ((PRECISION) 1)
 #undef TWO
-#define TWO ((ZOLO_PRECISION) 2)
+#define TWO ((PRECISION) 2)
 
 /* Evaluate the rational approximation R(x) using the factored form */
 
-static ZOLO_PRECISION zolotarev_eval(ZOLO_PRECISION x, zolotarev_data* rdata) {
+static PRECISION zolotarev_eval(PRECISION x, zolotarev_data* rdata) {
   int m;
-  ZOLO_PRECISION R;
+  PRECISION R;
 
   if (rdata -> type == 0) {
     R = rdata -> A * x;
@@ -551,9 +551,9 @@ static ZOLO_PRECISION zolotarev_eval(ZOLO_PRECISION x, zolotarev_data* rdata) {
 
 /* Evaluate the rational approximation R(x) using the partial fraction form */
 
-static ZOLO_PRECISION zolotarev_partfrac_eval(ZOLO_PRECISION x, zolotarev_data* rdata) {
+static PRECISION zolotarev_partfrac_eval(PRECISION x, zolotarev_data* rdata) {
   int m;
-  ZOLO_PRECISION R = rdata -> alpha[rdata -> da - 1];
+  PRECISION R = rdata -> alpha[rdata -> da - 1];
   for (m = 0; m < rdata -> dd; m++)
     R += rdata -> alpha[m] / (x * x - rdata -> ap[m]);
   if (rdata -> type == 1) R += rdata -> alpha[rdata -> dd] / (x * x);
@@ -568,18 +568,18 @@ static ZOLO_PRECISION zolotarev_partfrac_eval(ZOLO_PRECISION x, zolotarev_data*
  * non-signalling overflow this will work correctly since 1/(1/0) = 1/INF = 0,
  * but with signalling overflow you will get an error message. */
 
-static ZOLO_PRECISION zolotarev_contfrac_eval(ZOLO_PRECISION x, zolotarev_data* rdata) {
+static PRECISION zolotarev_contfrac_eval(PRECISION x, zolotarev_data* rdata) {
   int m;
-  ZOLO_PRECISION R = rdata -> beta[0] * x;
+  PRECISION R = rdata -> beta[0] * x;
   for (m = 1; m < rdata -> db; m++) R = rdata -> beta[m] * x + ONE / R;
   return R;
 }    
 
 /* Evaluate the rational approximation R(x) using Cayley form */
 
-static ZOLO_PRECISION zolotarev_cayley_eval(ZOLO_PRECISION x, zolotarev_data* rdata) {
+static PRECISION zolotarev_cayley_eval(PRECISION x, zolotarev_data* rdata) {
   int m;
-  ZOLO_PRECISION T;
+  PRECISION T;
 
   T = rdata -> type == 0 ? ONE : -ONE;
   for (m = 0; m < rdata -> n; m++)
@@ -607,7 +607,7 @@ int main(int argc, char** argv) {
   int m, n, plotpts = 5000, type = 0;
   float eps, x, ypferr, ycferr, ycaylerr, maxypferr, maxycferr, maxycaylerr;
   zolotarev_data *rdata;
-  ZOLO_PRECISION y;
+  PRECISION y;
   FILE *plot_function, *plot_error, 
     *plot_partfrac, *plot_contfrac, *plot_cayley;
 
@@ -626,13 +626,13 @@ int main(int argc, char** argv) {
   }
 
   rdata = type == 2 
-    ? higham((ZOLO_PRECISION) eps, n) 
-    : zolotarev((ZOLO_PRECISION) eps, n, type);
+    ? higham((PRECISION) eps, n) 
+    : zolotarev((PRECISION) eps, n, type);
 
   printf("Zolotarev Test: R(epsilon = %g, n = %d, type = %d)\n\t" 
 	 STRINGIFY(VERSION) "\n\t" STRINGIFY(HVERSION)
 	 "\n\tINTERNAL_PRECISION = " STRINGIFY(INTERNAL_PRECISION)
-	 "\tZOLO_PRECISION = " STRINGIFY(ZOLO_PRECISION)
+	 "\tPRECISION = " STRINGIFY(PRECISION)
 	 "\n\n\tRational approximation of degree (%d,%d), %s at x = 0\n"
 	 "\tDelta = %g (maximum error)\n\n"
 	 "\tA = %g (overall factor)\n",
@@ -681,15 +681,15 @@ int main(int argc, char** argv) {
     x = 2.4 * (float) m / plotpts - 1.2;
     if (rdata -> type == 0 || fabs(x) * (float) plotpts > 1.0) {
       /* skip x = 0 for type 1, as R(0) is singular */
-      y = zolotarev_eval((ZOLO_PRECISION) x, rdata);
+      y = zolotarev_eval((PRECISION) x, rdata);
       fprintf(plot_function, "%g %g\n", x, (float) y);
       fprintf(plot_error, "%g %g\n",
 	      x, (float)((y - ((x > 0.0 ? ONE : -ONE))) / rdata -> Delta));
-      ypferr = (float)((zolotarev_partfrac_eval((ZOLO_PRECISION) x, rdata) - y)
+      ypferr = (float)((zolotarev_partfrac_eval((PRECISION) x, rdata) - y)
 		       / rdata -> Delta);
-      ycferr = (float)((zolotarev_contfrac_eval((ZOLO_PRECISION) x, rdata) - y)
+      ycferr = (float)((zolotarev_contfrac_eval((PRECISION) x, rdata) - y)
 		       / rdata -> Delta);
-      ycaylerr = (float)((zolotarev_cayley_eval((ZOLO_PRECISION) x, rdata) - y)
+      ycaylerr = (float)((zolotarev_cayley_eval((PRECISION) x, rdata) - y)
 		       / rdata -> Delta);
       if (fabs(x) < 1.0 && fabs(x) > rdata -> epsilon) {
 	maxypferr = MAX(maxypferr, fabs(ypferr));

Grid/algorithms/approx/Zolotarev.h

@@ -9,10 +9,10 @@ NAMESPACE_BEGIN(Approx);
 #define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
 
 #ifndef ZOLOTAREV_INTERNAL
-#ifndef ZOLO_PRECISION
-#define ZOLO_PRECISION double
+#ifndef PRECISION
+#define PRECISION double
 #endif
-#define ZPRECISION ZOLO_PRECISION
+#define ZPRECISION PRECISION
 #define ZOLOTAREV_DATA zolotarev_data
 #endif
 
@@ -77,8 +77,8 @@ typedef struct {
  * zolotarev_data structure. The arguments must satisfy the constraints that
  * epsilon > 0, n > 0, and type = 0 or 1. */
 
-ZOLOTAREV_DATA* higham(ZOLO_PRECISION epsilon, int n) ;
-ZOLOTAREV_DATA* zolotarev(ZOLO_PRECISION epsilon, int n, int type);
+ZOLOTAREV_DATA* higham(PRECISION epsilon, int n) ;
+ZOLOTAREV_DATA* zolotarev(PRECISION epsilon, int n, int type);
 void zolotarev_free(zolotarev_data *zdata);
 #endif
 
@@ -86,4 +86,3 @@ void zolotarev_free(zolotarev_data *zdata);
 NAMESPACE_END(Approx);
 NAMESPACE_END(Grid);
 #endif
-

Grid/algorithms/blas/BatchedBlas.cc

@@ -1,34 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: BatchedBlas.h
-
-    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/GridCore.h>
-#include <Grid/algorithms/blas/BatchedBlas.h>
-NAMESPACE_BEGIN(Grid);
-gridblasHandle_t GridBLAS::gridblasHandle;
-int              GridBLAS::gridblasInit;
-NAMESPACE_END(Grid);
-

Grid/algorithms/blas/BatchedBlas.h

@@ -1,727 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: BatchedBlas.h
-
-    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 */
-#pragma once
-
-#ifdef GRID_HIP
-#include <hipblas/hipblas.h>
-#endif
-#ifdef GRID_CUDA
-#include <cublas_v2.h>
-#endif
-#ifdef GRID_SYCL
-#include <oneapi/mkl.hpp>
-#endif
-#if 0
-#define GRID_ONE_MKL
-#endif
-#ifdef GRID_ONE_MKL
-#include <oneapi/mkl.hpp>
-#endif
-///////////////////////////////////////////////////////////////////////	  
-// Need to rearrange lattice data to be in the right format for a
-// batched multiply. Might as well make these static, dense packed
-///////////////////////////////////////////////////////////////////////
-NAMESPACE_BEGIN(Grid);
-#ifdef GRID_HIP
-  typedef hipblasHandle_t gridblasHandle_t;
-#endif
-#ifdef GRID_CUDA
-  typedef cublasHandle_t gridblasHandle_t;
-#endif
-#ifdef GRID_SYCL
-  typedef cl::sycl::queue *gridblasHandle_t;
-#endif
-#ifdef GRID_ONE_MKL
-  typedef cl::sycl::queue *gridblasHandle_t;
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP) && !defined(GRID_ONE_MKL)
-  typedef int32_t gridblasHandle_t;
-#endif
-
-enum GridBLASOperation_t { GridBLAS_OP_N, GridBLAS_OP_T, GridBLAS_OP_C } ;
-
-class GridBLAS {
-public:
-
-  
-  static gridblasHandle_t gridblasHandle;
-  static int            gridblasInit;
-  
-  static void Init(void)
-  {
-    if ( ! gridblasInit ) {
-#ifdef GRID_CUDA
-      std::cout << "cublasCreate"<<std::endl;
-      cublasCreate(&gridblasHandle);
-      cublasSetPointerMode(gridblasHandle, CUBLAS_POINTER_MODE_DEVICE);
-#endif
-#ifdef GRID_HIP
-      std::cout << "hipblasCreate"<<std::endl;
-      hipblasCreate(&gridblasHandle);
-#endif
-#ifdef GRID_SYCL
-      gridblasHandle = theGridAccelerator;
-#endif
-#ifdef GRID_ONE_MKL
-      cl::sycl::cpu_selector selector;
-      cl::sycl::device selectedDevice { selector };
-      gridblasHandle =new sycl::queue (selectedDevice);
-#endif
-      gridblasInit=1;
-    }
-  }
-  
-  // Force construct once
-  GridBLAS() { Init(); };
-  ~GridBLAS() { };
-  
-  /////////////////////////////////////////////////////////////////////////////////////
-  // BLAS GEMM conventions:
-  /////////////////////////////////////////////////////////////////////////////////////
-  // - C = alpha A * B + beta C
-  // Dimensions:
-  // - C_m.n
-  // - A_m.k
-  // - B_k.n
-  // - Flops = 8 M N K
-  // - Bytes = 2*sizeof(word) * (MN+MK+KN)
-  // M=60, N=12
-  // Flop/Byte = 8 . 60.60.12 / (60.12+60.60+60.12)/16 = 4 so expect about 4 TF/s on a GCD
-  /////////////////////////////////////////////////////////////////////////////////////
-  void synchronise(void)
-  {
-#ifdef GRID_HIP
-    auto err = hipDeviceSynchronize();
-    assert(err==hipSuccess);
-#endif
-#ifdef GRID_CUDA
-    auto err = cudaDeviceSynchronize();
-    assert(err==cudaSuccess);
-#endif
-#ifdef GRID_SYCL
-    accelerator_barrier();
-#endif
-#ifdef GRID_ONE_MKL
-    gridblasHandle->wait();
-#endif
-  }
-  
-  void gemmBatched(int m,int n, int k,
-		   ComplexD alpha,
-		   deviceVector<ComplexD*> &Amk,  // pointer list to matrices
-		   deviceVector<ComplexD*> &Bkn,
-		   ComplexD beta,
-		   deviceVector<ComplexD*> &Cmn)
-  {
-    gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
-		m,n,k,
-		alpha,
-		Amk,
-		Bkn,
-		beta,
-		Cmn);
-  }
-  void gemmBatched(int m,int n, int k,
-		   ComplexF alpha,
-		   deviceVector<ComplexF*> &Amk,  // pointer list to matrices
-		   deviceVector<ComplexF*> &Bkn,
-		   ComplexF beta,
-		   deviceVector<ComplexF*> &Cmn)
-  {
-    gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
-		m,n,k,
-		alpha,
-		Amk,
-		Bkn,
-		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,
-		   int m,int n, int k,
-		   ComplexD alpha,
-		   deviceVector<ComplexD*> &Amk,  // pointer list to matrices
-		   deviceVector<ComplexD*> &Bkn,
-		   ComplexD beta,
-		   deviceVector<ComplexD*> &Cmn)
-  {
-    RealD t2=usecond();
-    int32_t batchCount = Amk.size();
-    assert(Bkn.size()==batchCount);
-    assert(Cmn.size()==batchCount);
-
-    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();
-    //    std::cout << "ZgemmBatched mnk  "<<m<<","<<n<<","<<k<<" count "<<batchCount<<std::endl;
-#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 = hipblasZgemmBatched(gridblasHandle,
-				   hOpA,
-				   hOpB,
-				   m,n,k,
-				   (hipblasDoubleComplex *) &alpha_p[0],
-				   (hipblasDoubleComplex **)&Amk[0], lda,
-				   (hipblasDoubleComplex **)&Bkn[0], ldb,
-				   (hipblasDoubleComplex *) &beta_p[0],
-				   (hipblasDoubleComplex **)&Cmn[0], ldc,
-				   batchCount);
-    //	 std::cout << " hipblas return code " <<(int)err<<std::endl;
-    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 = cublasZgemmBatched(gridblasHandle,
-				  hOpA,
-				  hOpB,
-				  m,n,k,
-				  (cuDoubleComplex *) &alpha_p[0],
-				  (cuDoubleComplex **)&Amk[0], lda,
-				  (cuDoubleComplex **)&Bkn[0], ldb,
-				  (cuDoubleComplex *) &beta_p[0],
-				  (cuDoubleComplex **)&Cmn[0], ldc,
-				  batchCount);
-    assert(err==CUBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_SYCL
-    //MKL’s cblas_<T>gemm_batch & OneAPI
-#warning "oneMKL implementation not built "
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    // Need a default/reference implementation
-    int sda = lda*k;
-    int sdb = ldb*k;
-    int sdc = ldc*n;
-    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[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	}
-      }
-    }
-#endif
-    //    synchronise();
-     RealD t1=usecond();
-     RealD flops = 8.0*m*n*k*batchCount;
-     RealD bytes = 1.0*sizeof(ComplexD)*(m*k+k*n+m*n)*batchCount;
-     //     std::cout <<GridLogMessage<< " batched Blas copy "<<(t0-t2)/1.e3 <<" ms "<<std::endl;
-     //     std::cout <<GridLogMessage<< " batched Blas zGemm call "<<m<<","<<n<<","<<k<<" "<< flops/(t1-t0)/1.e3 <<" GF/s "<<(t1-t0)/1.e3<<" ms "<<std::endl;
-     //     std::cout <<GridLogMessage<< " batched Blas zGemm call "<<m<<","<<n<<","<<k<<" "<< bytes/(t1-t0)/1.e3 <<" GB/s "<<(t1-t0)/1.e3<<" ms "<<std::endl;
-  }
-
-  void gemmBatched(GridBLASOperation_t OpA,
-		   GridBLASOperation_t OpB,
-		   int m,int n, int k,
-		   ComplexF alpha,
-		   deviceVector<ComplexF*> &Amk,  // pointer list to matrices
-		   deviceVector<ComplexF*> &Bkn,
-		   ComplexF beta,
-		   deviceVector<ComplexF*> &Cmn)
-  {
-    RealD t2=usecond();
-    int32_t batchCount = Amk.size();
-
-    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();
-
-    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 = hipblasCgemmBatched(gridblasHandle,
-				   hOpA,
-				   hOpB,
-				   m,n,k,
-				   (hipblasComplex *) &alpha_p[0],
-				   (hipblasComplex **)&Amk[0], lda,
-				   (hipblasComplex **)&Bkn[0], ldb,
-				   (hipblasComplex *) &beta_p[0],
-				   (hipblasComplex **)&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 = cublasCgemmBatched(gridblasHandle,
-				  hOpA,
-				  hOpB,
-				  m,n,k,
-				  (cuComplex *) &alpha_p[0],
-				  (cuComplex **)&Amk[0], lda,
-				  (cuComplex **)&Bkn[0], ldb,
-				  (cuComplex *) &beta_p[0],
-				  (cuComplex **)&Cmn[0], ldc,
-				  batchCount);
-    assert(err==CUBLAS_STATUS_SUCCESS);
-#endif
-#ifdef GRID_SYCL
-    //MKL’s cblas_<T>gemm_batch & OneAPI
-#warning "oneMKL implementation not built "
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    int sda = lda*k;
-    int sdb = ldb*k;
-    int sdc = ldc*n;
-    ComplexF alphaf(real(alpha),imag(alpha));
-    ComplexF betaf(real(beta),imag(beta));
-    // 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) {
-	  ComplexF c_mn(0.0);
-	  for (int kk = 0; kk < k; ++kk)
-	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Cmn[p][mm + nn*ldc] =  (alphaf)*c_mn + (betaf)*Cmn[p][mm + nn*ldc ];
-	}
-      }
-    }
-#endif
-     RealD t1=usecond();
-     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();
-
-    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
-    //MKL’s cblas_<T>gemm_batch & OneAPI
-#warning "oneMKL implementation not built "
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    int sda = lda*k;
-    int sdb = ldb*k;
-    int sdc = ldc*n;
-    // 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) {
-	  RealD c_mn(0.0);
-	  for (int kk = 0; kk < k; ++kk)
-	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	}
-      }
-    }
-#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();
-
-    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 batchCount64=batchCount;
-      oneapi::mkl::blas::column_major::gemm_batch(*theGridAccelerator,
-      onemkl::transpose::N,
-      onemkl::transpose::N,
-      &m64,&n64,&k64,
-      (double *) &alpha_p[0],
-      (double **)&Amk[0], lda,
-      (double **)&Bkn[0], ldb,
-      (double *) &beta_p[0],
-      (double **)&Cmn[0], ldc,
-      1,&batchCount64);
-     */
-    //MKL’s cblas_<T>gemm_batch & OneAPI
-#warning "oneMKL implementation not built "
-#endif
-#if !defined(GRID_SYCL) && !defined(GRID_CUDA) && !defined(GRID_HIP)
-    int sda = lda*k;
-    int sdb = ldb*k;
-    int sdc = ldc*n;
-    // 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) {
-	  RealD c_mn(0.0);
-	  for (int kk = 0; kk < k; ++kk)
-	    c_mn += Amk[p][mm + kk*lda ] * Bkn[p][kk + nn*ldb];
-	  Cmn[p][mm + nn*ldc] =  (alpha)*c_mn + (beta)*Cmn[p][mm + nn*ldc ];
-	}
-      }
-    }
-#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);
-#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
-  }
-
-  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<ComplexD> A(N_A); acceleratorMemSet(&A[0],0,N_A*sizeof(ComplexD));
-    deviceVector<ComplexD> B(N_B); acceleratorMemSet(&B[0],0,N_B*sizeof(ComplexD));
-    deviceVector<ComplexD> C(N_C); acceleratorMemSet(&C[0],0,N_C*sizeof(ComplexD));
-    ComplexD alpha(1.0);
-    ComplexD beta (1.0);
-    RealD flops = 8.0*M*N*K*BATCH;
-    int ncall=10;
-    RealD t0 = usecond();
-    for(int i=0;i<ncall;i++){
-      gemmStridedBatched(M,N,K,
-			 alpha,
-			 &A[0], // m x k 
-			 &B[0], // k x n
-			 beta, 
-			 &C[0], // m x n
-			 BATCH);
-    }
-    synchronise();
-    RealD t1 = usecond();
-    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/allocator/AlignedAllocator.h

@@ -176,7 +176,6 @@ template<class T> using cshiftAllocator = std::allocator<T>;
 template<class T> using Vector        = std::vector<T,uvmAllocator<T> >;           
 template<class T> using stencilVector = std::vector<T,alignedAllocator<T> >;           
 template<class T> using commVector = std::vector<T,devAllocator<T> >;
-template<class T> using deviceVector  = std::vector<T,devAllocator<T> >;
 template<class T> using cshiftVector = std::vector<T,cshiftAllocator<T> >;
 
 NAMESPACE_END(Grid);

Grid/communicator/Communicator_mpi3.cc

@@ -348,7 +348,6 @@ double CartesianCommunicator::StencilSendToRecvFrom( void *xmit,
   return offbytes;
 }
 
-#undef NVLINK_GET // Define to use get instead of put DMA
 double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsRequest_t> &list,
 							 void *xmit,
 							 int dest,int dox,
@@ -381,15 +380,9 @@ double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsReques
       list.push_back(rrq);
       off_node_bytes+=rbytes;
     }
-#ifdef NVLINK_GET
-      void *shm = (void *) this->ShmBufferTranslate(from,xmit);
-      assert(shm!=NULL);
-      acceleratorCopyDeviceToDeviceAsynch(shm,recv,rbytes);
-#endif
   }
   
   if (dox) {
-    //  rcrc = crc32(rcrc,(unsigned char *)recv,bytes);
     if ( (gdest == MPI_UNDEFINED) || Stencil_force_mpi ) {
       tag= dir+_processor*32;
       ierr =MPI_Isend(xmit, xbytes, MPI_CHAR,dest,tag,communicator_halo[commdir],&xrq);
@@ -397,12 +390,9 @@ double CartesianCommunicator::StencilSendToRecvFromBegin(std::vector<CommsReques
       list.push_back(xrq);
       off_node_bytes+=xbytes;
     } else {
-#ifndef NVLINK_GET
       void *shm = (void *) this->ShmBufferTranslate(dest,recv);
       assert(shm!=NULL);
       acceleratorCopyDeviceToDeviceAsynch(xmit,shm,xbytes);
-#endif
-      
     }
   }
 
@@ -412,8 +402,6 @@ void CartesianCommunicator::StencilSendToRecvFromComplete(std::vector<CommsReque
 {
   int nreq=list.size();
 
-  acceleratorCopySynchronise();
-
   if (nreq==0) return;
 
   std::vector<MPI_Status> status(nreq);

Grid/communicator/SharedMemory.cc

@@ -40,9 +40,6 @@ int                 GlobalSharedMemory::_ShmAlloc;
 uint64_t            GlobalSharedMemory::_ShmAllocBytes;
 
 std::vector<void *> GlobalSharedMemory::WorldShmCommBufs;
-#ifndef ACCELERATOR_AWARE_MPI
-void * GlobalSharedMemory::HostCommBuf;
-#endif
 
 Grid_MPI_Comm       GlobalSharedMemory::WorldShmComm;
 int                 GlobalSharedMemory::WorldShmRank;
@@ -69,26 +66,6 @@ void GlobalSharedMemory::SharedMemoryFree(void)
 /////////////////////////////////
 // Alloc, free shmem region
 /////////////////////////////////
-#ifndef ACCELERATOR_AWARE_MPI
-void *SharedMemory::HostBufferMalloc(size_t bytes){
-  void *ptr = (void *)host_heap_top;
-  host_heap_top  += bytes;
-  host_heap_bytes+= bytes;
-  if (host_heap_bytes >= host_heap_size) {
-    std::cout<< " HostBufferMalloc exceeded heap size -- try increasing with --shm <MB> flag" <<std::endl;
-    std::cout<< " Parameter specified in units of MB (megabytes) " <<std::endl;
-    std::cout<< " Current alloc is " << (bytes/(1024*1024)) <<"MB"<<std::endl;
-    std::cout<< " Current bytes is " << (host_heap_bytes/(1024*1024)) <<"MB"<<std::endl;
-    std::cout<< " Current heap  is " << (host_heap_size/(1024*1024)) <<"MB"<<std::endl;
-    assert(host_heap_bytes<host_heap_size);
-  }
-  return ptr;
-}
-void SharedMemory::HostBufferFreeAll(void) { 
-  host_heap_top  =(size_t)HostCommBuf;
-  host_heap_bytes=0;
-}
-#endif
 void *SharedMemory::ShmBufferMalloc(size_t bytes){
   //  bytes = (bytes+sizeof(vRealD))&(~(sizeof(vRealD)-1));// align up bytes
   void *ptr = (void *)heap_top;

Grid/communicator/SharedMemory.h

@@ -75,9 +75,7 @@ public:
   static int           Hugepages;
 
   static std::vector<void *> WorldShmCommBufs;
-#ifndef ACCELERATOR_AWARE_MPI
-  static void *HostCommBuf;
-#endif
+
   static Grid_MPI_Comm WorldComm;
   static int           WorldRank;
   static int           WorldSize;
@@ -122,13 +120,6 @@ private:
   size_t heap_bytes;
   size_t heap_size;
 
-#ifndef ACCELERATOR_AWARE_MPI
-  size_t host_heap_top;  // set in free all
-  size_t host_heap_bytes;// set in free all
-  void *HostCommBuf;     // set in SetCommunicator
-  size_t host_heap_size; // set in SetCommunicator
-#endif
-  
 protected:
 
   Grid_MPI_Comm    ShmComm; // for barriers
@@ -160,10 +151,7 @@ public:
   void *ShmBufferTranslate(int rank,void * local_p);
   void *ShmBufferMalloc(size_t bytes);
   void  ShmBufferFreeAll(void) ;
-#ifndef ACCELERATOR_AWARE_MPI
-  void *HostBufferMalloc(size_t bytes);
-  void HostBufferFreeAll(void);
-#endif  
+  
   //////////////////////////////////////////////////////////////////////////
   // Make info on Nodes & ranks and Shared memory available
   //////////////////////////////////////////////////////////////////////////

Grid/communicator/SharedMemoryMPI.cc

@@ -39,11 +39,9 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #include <hip/hip_runtime_api.h>
 #endif
 #ifdef GRID_SYCL
-#ifdef ACCELERATOR_AWARE_MPI
 #define GRID_SYCL_LEVEL_ZERO_IPC
-#define SHM_SOCKETS
-#endif 
 #include <syscall.h>
+#define SHM_SOCKETS 
 #endif
 
 #include <sys/socket.h>
@@ -514,6 +512,46 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 // Hugetlbfs mapping intended
 ////////////////////////////////////////////////////////////////////////////////////////////
 #if defined(GRID_CUDA) ||defined(GRID_HIP)  || defined(GRID_SYCL)
+
+//if defined(GRID_SYCL)
+#if 0
+void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
+{
+  void * ShmCommBuf ; 
+  assert(_ShmSetup==1);
+  assert(_ShmAlloc==0);
+
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // allocate the pointer array for shared windows for our group
+  //////////////////////////////////////////////////////////////////////////////////////////////////////////
+  MPI_Barrier(WorldShmComm);
+  WorldShmCommBufs.resize(WorldShmSize);
+
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
+  // Each MPI rank should allocate our own buffer
+  ///////////////////////////////////////////////////////////////////////////////////////////////////////////
+  ShmCommBuf = acceleratorAllocDevice(bytes);
+
+  if (ShmCommBuf == (void *)NULL ) {
+    std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
+    exit(EXIT_FAILURE);  
+  }
+
+  std::cout << WorldRank << Mheader " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
+	    << "bytes at "<< std::hex<< ShmCommBuf <<std::dec<<" for comms buffers " <<std::endl;
+
+  SharedMemoryZero(ShmCommBuf,bytes);
+
+  assert(WorldShmSize == 1);
+  for(int r=0;r<WorldShmSize;r++){
+    WorldShmCommBufs[r] = ShmCommBuf;
+  }
+  _ShmAllocBytes=bytes;
+  _ShmAlloc=1;
+}
+#endif
+
+#if defined(GRID_CUDA) ||defined(GRID_HIP) ||defined(GRID_SYCL)  
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
   void * ShmCommBuf ; 
@@ -536,9 +574,6 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
   // Each MPI rank should allocate our own buffer
   ///////////////////////////////////////////////////////////////////////////////////////////////////////////
-#ifndef ACCELERATOR_AWARE_MPI
-  HostCommBuf= malloc(bytes);
-#endif  
   ShmCommBuf = acceleratorAllocDevice(bytes);
   if (ShmCommBuf == (void *)NULL ) {
     std::cerr << " SharedMemoryMPI.cc acceleratorAllocDevice failed NULL pointer for " << bytes<<" bytes " << std::endl;
@@ -703,6 +738,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
   _ShmAllocBytes=bytes;
   _ShmAlloc=1;
 }
+#endif
 
 #else 
 #ifdef GRID_MPI3_SHMMMAP
@@ -926,12 +962,6 @@ void SharedMemory::SetCommunicator(Grid_MPI_Comm comm)
   }
   ShmBufferFreeAll();
 
-#ifndef ACCELERATOR_AWARE_MPI
-  host_heap_size = heap_size;
-  HostCommBuf= GlobalSharedMemory::HostCommBuf;
-  HostBufferFreeAll();
-#endif  
-
   /////////////////////////////////////////////////////////////////////
   // find comm ranks in our SHM group (i.e. which ranks are on our node)
   /////////////////////////////////////////////////////////////////////

Grid/lattice/Lattice.h

@@ -35,7 +35,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_transpose.h>
 #include <Grid/lattice/Lattice_local.h>
 #include <Grid/lattice/Lattice_reduction.h>
-#include <Grid/lattice/Lattice_crc.h>
 #include <Grid/lattice/Lattice_peekpoke.h>
 #include <Grid/lattice/Lattice_reality.h>
 #include <Grid/lattice/Lattice_real_imag.h>
@@ -47,4 +46,5 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 #include <Grid/lattice/Lattice_unary.h>
 #include <Grid/lattice/Lattice_transfer.h>
 #include <Grid/lattice/Lattice_basis.h>
+#include <Grid/lattice/Lattice_crc.h>
 #include <Grid/lattice/PaddedCell.h>

Grid/lattice/Lattice_crc.h

@@ -42,13 +42,13 @@ template<class vobj> void DumpSliceNorm(std::string s,Lattice<vobj> &f,int mu=-1
   }
 }
 
-template<class vobj> uint32_t crc(const Lattice<vobj> & buf)
+template<class vobj> uint32_t crc(Lattice<vobj> & buf)
 {
   autoView( buf_v , buf, CpuRead);
   return ::crc32(0L,(unsigned char *)&buf_v[0],(size_t)sizeof(vobj)*buf.oSites());
 }
 
-#define CRC(U) std::cerr << "FingerPrint "<<__FILE__ <<" "<< __LINE__ <<" "<< #U <<" "<<crc(U)<<std::endl;
+#define CRC(U) std::cout << "FingerPrint "<<__FILE__ <<" "<< __LINE__ <<" "<< #U <<" "<<crc(U)<<std::endl;
 
 NAMESPACE_END(Grid);
 

Grid/lattice/Lattice_reduction.h

@@ -281,29 +281,11 @@ inline ComplexD rankInnerProduct(const Lattice<vobj> &left,const Lattice<vobj> &
   return nrm;
 }
 
-
 template<class vobj>
 inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) {
   GridBase *grid = left.Grid();
-
-#ifdef GRID_SYCL
-  uint64_t csum=0;
-  if ( FlightRecorder::LoggingMode != FlightRecorder::LoggingModeNone)
-  {
-    // Hack
-    // Fast integer xor checksum. Can also be used in comms now.
-    autoView(l_v,left,AcceleratorRead);
-    Integer words = left.Grid()->oSites()*sizeof(vobj)/sizeof(uint64_t);
-    uint64_t *base= (uint64_t *)&l_v[0];
-    csum=svm_xor(base,words);
-  }
-  FlightRecorder::CsumLog(csum);
-#endif
   ComplexD nrm = rankInnerProduct(left,right);
-  RealD local = real(nrm);
-  FlightRecorder::NormLog(real(nrm)); 
   grid->GlobalSum(nrm);
-  FlightRecorder::ReductionLog(local,real(nrm)); 
   return nrm;
 }
 

Grid/lattice/Lattice_reduction_sycl.h

@@ -69,29 +69,28 @@ inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osite
   return result;
 }
 
-
-template<class Word> Word svm_xor(Word *vec,uint64_t L)
-{
-  Word xorResult; xorResult = 0;
-  Word *d_sum =(Word *)cl::sycl::malloc_shared(sizeof(Word),*theGridAccelerator);
-  Word identity;  identity=0;
-  theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
-     auto Reduction = cl::sycl::reduction(d_sum,identity,std::bit_xor<>());
-     cgh.parallel_for(cl::sycl::range<1>{L},
-		      Reduction,
-		      [=] (cl::sycl::id<1> index, auto &sum) {
-	 sum ^=vec[index];
-     });
-   });
-  theGridAccelerator->wait();
-  Word ret = d_sum[0];
-  free(d_sum,*theGridAccelerator);
-  return ret;
-}
-
 NAMESPACE_END(Grid);
 
 /*
+template<class Double> Double svm_reduce(Double *vec,uint64_t L)
+{
+  Double sumResult; zeroit(sumResult);
+  Double *d_sum =(Double *)cl::sycl::malloc_shared(sizeof(Double),*theGridAccelerator);
+  Double identity;  zeroit(identity);
+  theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
+     auto Reduction = cl::sycl::reduction(d_sum,identity,std::plus<>());
+     cgh.parallel_for(cl::sycl::range<1>{L},
+		      Reduction,
+		      [=] (cl::sycl::id<1> index, auto &sum) {
+	 sum +=vec[index];
+     });
+   });
+  theGridAccelerator->wait();
+  Double ret = d_sum[0];
+  free(d_sum,*theGridAccelerator);
+  std::cout << " svm_reduce finished "<<L<<" sites sum = " << ret <<std::endl;
+  return ret;
+}
 
 template <class vobj>
 inline typename vobj::scalar_objectD sumD_gpu_repack(const vobj *lat, Integer osites)

Grid/lattice/Lattice_rng.h

@@ -411,7 +411,7 @@ public:
       std::cout << GridLogMessage << "Seed SHA256: " << GridChecksum::sha256_string(seeds) << std::endl;
       SeedFixedIntegers(seeds);
     }
-  void SeedFixedIntegers(const std::vector<int> &seeds, int britney=0){
+  void SeedFixedIntegers(const std::vector<int> &seeds){
 
     // Everyone generates the same seed_seq based on input seeds
     CartesianCommunicator::BroadcastWorld(0,(void *)&seeds[0],sizeof(int)*seeds.size());
@@ -428,6 +428,7 @@ public:
     // MT implementation does not implement fast discard even though
     // in principle this is possible
     ////////////////////////////////////////////////
+#if 1
     thread_for( lidx, _grid->lSites(), {
 
 	int gidx;
@@ -448,12 +449,29 @@ public:
 	
 	int l_idx=generator_idx(o_idx,i_idx);
 	_generators[l_idx] = master_engine;
-	if ( britney ) { 
-	  Skip(_generators[l_idx],l_idx); // Skip to next RNG sequence
-	} else { 	
+	Skip(_generators[l_idx],gidx); // Skip to next RNG sequence
+    });
+#else
+    // Everybody loops over global volume.
+    thread_for( gidx, _grid->_gsites, {
+
+	// Where is it?
+	int rank;
+	int o_idx;
+	int i_idx;
+
+	Coordinate gcoor;
+	_grid->GlobalIndexToGlobalCoor(gidx,gcoor);
+	_grid->GlobalCoorToRankIndex(rank,o_idx,i_idx,gcoor);
+	
+	// If this is one of mine we take it
+	if( rank == _grid->ThisRank() ){
+	  int l_idx=generator_idx(o_idx,i_idx);
+	  _generators[l_idx] = master_engine;
 	  Skip(_generators[l_idx],gidx); // Skip to next RNG sequence
 	}
     });
+#endif
 #else 
     ////////////////////////////////////////////////////////////////
     // Machine and thread decomposition dependent seeding is efficient

Grid/qcd/action/fermion/implementation/StaggeredKernelsImplementation.h

@@ -280,16 +280,20 @@ void StaggeredKernels<Impl>::DhopImproved(StencilImpl &st, LebesgueOrder &lo,
 
   if( interior && exterior ) { 
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGeneric,1); return;}
-    if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHand,1);    return;}
 #ifndef GRID_CUDA
+    if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHand,1);    return;}
     if (Opt == OptInlineAsm  ) {  ASM_CALL(DhopSiteAsm);     return;}
 #endif
   } else if( interior ) {
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGenericInt,1); return;}
+#ifndef GRID_CUDA
     if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHandInt,1);    return;}
+#endif
   } else if( exterior ) { 
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGenericExt,1); return;}
+#ifndef GRID_CUDA
     if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHandExt,1);    return;}
+#endif
   }
   assert(0 && " Kernel optimisation case not covered ");
 }
@@ -318,13 +322,19 @@ void StaggeredKernels<Impl>::DhopNaive(StencilImpl &st, LebesgueOrder &lo,
   
   if( interior && exterior ) { 
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGeneric,0); return;}
+#ifndef GRID_CUDA
     if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHand,0);    return;}
+#endif
   } else if( interior ) {
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGenericInt,0); return;}
+#ifndef GRID_CUDA
     if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHandInt,0);    return;}
+#endif
   } else if( exterior ) { 
     if (Opt == OptGeneric    ) { KERNEL_CALL(DhopSiteGenericExt,0); return;}
+#ifndef GRID_CUDA
     if (Opt == OptHandUnroll ) { KERNEL_CALL(DhopSiteHandExt,0);    return;}
+#endif
   }
 }
 

Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h

@@ -462,7 +462,6 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
     autoView(st_v , st,AcceleratorRead);
 
    if( interior && exterior ) {
-     acceleratorFenceComputeStream();
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSite); return;}
      if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSite);    return;}
 #ifndef GRID_CUDA
@@ -496,7 +495,6 @@ void WilsonKernels<Impl>::DhopKernel(int Opt,StencilImpl &st,  DoubledGaugeField
     autoView(st_v ,st,AcceleratorRead);
 
    if( interior && exterior ) {
-     acceleratorFenceComputeStream();
      if (Opt == WilsonKernelsStatic::OptGeneric    ) { KERNEL_CALL(GenericDhopSiteDag); return;}
      if (Opt == WilsonKernelsStatic::OptHandUnroll ) { KERNEL_CALL(HandDhopSiteDag);    return;}
 #ifndef GRID_CUDA

Grid/qcd/smearing/HISQSmearing.h

@@ -138,7 +138,7 @@ public:
         for(int nu=0;nu<Nd;nu++) {
             appendShift(shifts,mu);
             appendShift(shifts,nu);
-            appendShift(shifts,shiftSignal::NO_SHIFT);
+            appendShift(shifts,NO_SHIFT);
             appendShift(shifts,mu,Back(nu));
             appendShift(shifts,Back(nu));
             appendShift(shifts,Back(mu));
@@ -173,6 +173,7 @@ public:
             auto gStencil_v = gStencil.View(); 
 
             accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 3-link constructs
+//            for(int site=0;site<Nsites;site++){ // ----------- 3-link constructs
                 stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
                 U3matrix U0, U1, U2, U3, U4, U5, W;
                 for(int nu=0;nu<Nd;nu++) {
@@ -215,6 +216,7 @@ public:
             })
 
             accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 5-link 
+//            for(int site=0;site<Nsites;site++){ // ----------- 5-link
                 stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
                 U3matrix U0, U1, U2, U3, U4, U5, W;
                 int sigmaIndex = 0;
@@ -252,6 +254,7 @@ public:
             })
 
             accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 7-link
+//            for(int site=0;site<Nsites;site++){ // ----------- 7-link
                 stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
                 U3matrix U0, U1, U2, U3, U4, U5, W;
                 int sigmaIndex = 0;

Grid/simd/Grid_vector_types.h

@@ -1133,13 +1133,4 @@ static_assert(sizeof(SIMD_Ftype) == sizeof(SIMD_Itype), "SIMD vector lengths inc
 
 NAMESPACE_END(Grid);
 
-#ifdef GRID_SYCL
-template<> struct sycl::is_device_copyable<Grid::vComplexF> : public std::true_type {};
-template<> struct sycl::is_device_copyable<Grid::vComplexD> : public std::true_type {};
-template<> struct sycl::is_device_copyable<Grid::vRealF   > : public std::true_type {};
-template<> struct sycl::is_device_copyable<Grid::vRealD   > : public std::true_type {};
-template<> struct sycl::is_device_copyable<Grid::vInteger > : public std::true_type {};
-#endif
-
-
 #endif

Grid/stencil/GeneralLocalStencil.h

@@ -141,14 +141,8 @@ public:
 ////////////////////////////////////////////////
 // Some machinery to streamline making a stencil 
 ////////////////////////////////////////////////
-
-class shiftSignal {
-public:
-    enum {
-        BACKWARD_CONST = 16,
-        NO_SHIFT       = -1
-    };
-};
+#define BACKWARD_CONST 16
+#define NO_SHIFT -1
 
 // TODO: put a check somewhere that BACKWARD_CONST > Nd!
 
@@ -156,16 +150,16 @@ public:
 inline int Back(const int dir) {
     // generalShift will use BACKWARD_CONST to determine whether we step forward or 
     // backward. Trick inspired by SIMULATeQCD. 
-    return dir + shiftSignal::BACKWARD_CONST;
+    return dir + BACKWARD_CONST;
 }
 
 /*!  @brief shift one unit in direction dir */
 template<typename... Args>
 void generalShift(Coordinate& shift, int dir) {
-    if (dir >= shiftSignal::BACKWARD_CONST) {
-        dir -= shiftSignal::BACKWARD_CONST;
+    if (dir >= BACKWARD_CONST) {
+        dir -= BACKWARD_CONST;
         shift[dir]+=-1;
-    } else if (dir == shiftSignal::NO_SHIFT) {
+    } else if (dir == NO_SHIFT) {
         ; // do nothing
     } else {
         shift[dir]+=1;
@@ -175,10 +169,10 @@ void generalShift(Coordinate& shift, int dir) {
 /*!  @brief follow a path of directions, shifting one unit in each direction */
 template<typename... Args>
 void generalShift(Coordinate& shift, int dir, Args... args) {
-    if (dir >= shiftSignal::BACKWARD_CONST) {
-        dir -= shiftSignal::BACKWARD_CONST;
+    if (dir >= BACKWARD_CONST) {
+        dir -= BACKWARD_CONST;
         shift[dir]+=-1;
-    } else if (dir == shiftSignal::NO_SHIFT) {
+    } else if (dir == NO_SHIFT) {
         ; // do nothing
     } else {
         shift[dir]+=1;

Grid/stencil/Stencil.h

@@ -70,6 +70,57 @@ struct DefaultImplParams {
 void Gather_plane_table_compute (GridBase *grid,int dimension,int plane,int cbmask,
 				 int off,std::vector<std::pair<int,int> > & table);
 
+/*
+template<class vobj,class cobj,class compressor>
+void Gather_plane_simple_table (commVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,cobj *buffer,compressor &compress, int off,int so)   __attribute__((noinline));
+
+template<class vobj,class cobj,class compressor>
+void Gather_plane_simple_table (commVector<std::pair<int,int> >& table,const Lattice<vobj> &rhs,cobj *buffer,compressor &compress, int off,int so)
+{
+  int num=table.size();
+  std::pair<int,int> *table_v = & table[0];
+
+  auto rhs_v = rhs.View(AcceleratorRead);
+  accelerator_forNB( i,num, vobj::Nsimd(), {
+    compress.Compress(buffer[off+table_v[i].first],rhs_v[so+table_v[i].second]);
+  });
+  rhs_v.ViewClose();
+}
+
+///////////////////////////////////////////////////////////////////
+// Gather for when there *is* need to SIMD split with compression
+///////////////////////////////////////////////////////////////////
+template<class cobj,class vobj,class compressor>
+void Gather_plane_exchange_table(const Lattice<vobj> &rhs,
+				 commVector<cobj *> pointers,
+				 int dimension,int plane,
+				 int cbmask,compressor &compress,int type) __attribute__((noinline));
+
+template<class cobj,class vobj,class compressor>
+void Gather_plane_exchange_table(commVector<std::pair<int,int> >& table,
+				 const Lattice<vobj> &rhs,
+				 std::vector<cobj *> &pointers,int dimension,int plane,int cbmask,
+				 compressor &compress,int type)
+{
+  assert( (table.size()&0x1)==0);
+  int num=table.size()/2;
+  int so  = plane*rhs.Grid()->_ostride[dimension]; // base offset for start of plane
+
+  auto rhs_v = rhs.View(AcceleratorRead);
+  auto rhs_p = &rhs_v[0];
+  auto p0=&pointers[0][0];
+  auto p1=&pointers[1][0];
+  auto tp=&table[0];
+  accelerator_forNB(j, num, vobj::Nsimd(), {
+      compress.CompressExchange(p0,p1, rhs_p, j,
+				so+tp[2*j  ].second,
+				so+tp[2*j+1].second,
+				type);
+  });
+  rhs_v.ViewClose();
+}
+*/
+
 void DslashResetCounts(void);
 void DslashGetCounts(uint64_t &dirichlet,uint64_t &partial,uint64_t &full);
 void DslashLogFull(void);
@@ -207,10 +258,6 @@ public:
   struct Packet {
     void * send_buf;
     void * recv_buf;
-#ifndef ACCELERATOR_AWARE_MPI
-    void * host_send_buf; // Allocate this if not MPI_CUDA_AWARE
-    void * host_recv_buf; // Allocate this if not MPI_CUDA_AWARE
-#endif
     Integer to_rank;
     Integer from_rank;
     Integer do_send;
@@ -277,7 +324,7 @@ public:
   Vector<int> surface_list;
 
   stencilVector<StencilEntry>  _entries; // Resident in managed memory
-  commVector<StencilEntry>     _entries_device; // Resident in device memory
+  commVector<StencilEntry>     _entries_device; // Resident in managed memory
   std::vector<Packet> Packets;
   std::vector<Merge> Mergers;
   std::vector<Merge> MergersSHM;
@@ -361,16 +408,33 @@ public:
   // Use OpenMP Tasks for cleaner ???
   // must be called *inside* parallel region
   //////////////////////////////////////////
+  /*
+  void CommunicateThreaded()
+  {
+#ifdef GRID_OMP
+    int mythread = omp_get_thread_num();
+    int nthreads = CartesianCommunicator::nCommThreads;
+#else
+    int mythread = 0;
+    int nthreads = 1;
+#endif
+    if (nthreads == -1) nthreads = 1;
+    if (mythread < nthreads) {
+      for (int i = mythread; i < Packets.size(); i += nthreads) {
+	uint64_t bytes = _grid->StencilSendToRecvFrom(Packets[i].send_buf,
+						      Packets[i].to_rank,
+						      Packets[i].recv_buf,
+						      Packets[i].from_rank,
+						      Packets[i].bytes,i);
+      }
+    }
+  }
+  */
   ////////////////////////////////////////////////////////////////////////
   // Non blocking send and receive. Necessarily parallel.
   ////////////////////////////////////////////////////////////////////////
   void CommunicateBegin(std::vector<std::vector<CommsRequest_t> > &reqs)
   {
-    // All GPU kernel tasks must complete
-    //    accelerator_barrier();     // All kernels should ALREADY be complete
-    //    _grid->StencilBarrier();   // Everyone is here, so noone running slow and still using receive buffer
-                               // But the HaloGather had a barrier too.
-#ifdef ACCELERATOR_AWARE_MPI
     for(int i=0;i<Packets.size();i++){
       _grid->StencilSendToRecvFromBegin(MpiReqs,
 					Packets[i].send_buf,
@@ -379,54 +443,16 @@ public:
 					Packets[i].from_rank,Packets[i].do_recv,
 					Packets[i].xbytes,Packets[i].rbytes,i);
     }
-#else
-#warning "Using COPY VIA HOST BUFFERS IN STENCIL"
-    for(int i=0;i<Packets.size();i++){
-      // Introduce a host buffer with a cheap slab allocator and zero cost wipe all
-      Packets[i].host_send_buf = _grid->HostBufferMalloc(Packets[i].xbytes);
-      Packets[i].host_recv_buf = _grid->HostBufferMalloc(Packets[i].rbytes);
-      if ( Packets[i].do_send ) {
-	acceleratorCopyFromDevice(Packets[i].send_buf, Packets[i].host_send_buf,Packets[i].xbytes);
-      }
-      _grid->StencilSendToRecvFromBegin(MpiReqs,
-					Packets[i].host_send_buf,
-					Packets[i].to_rank,Packets[i].do_send,
-					Packets[i].host_recv_buf,
-					Packets[i].from_rank,Packets[i].do_recv,
-					Packets[i].xbytes,Packets[i].rbytes,i);
-    }
-#endif
-    // Get comms started then run checksums
-    // Having this PRIOR to the dslash seems to make Sunspot work... (!)
-    for(int i=0;i<Packets.size();i++){
-      if ( Packets[i].do_send )
-	FlightRecorder::xmitLog(Packets[i].send_buf,Packets[i].xbytes);
-    }
   }
 
   void CommunicateComplete(std::vector<std::vector<CommsRequest_t> > &reqs)
   {
-    _grid->StencilSendToRecvFromComplete(MpiReqs,0); // MPI is done
+    _grid->StencilSendToRecvFromComplete(MpiReqs,0);
     if   ( this->partialDirichlet ) DslashLogPartial();
     else if ( this->fullDirichlet ) DslashLogDirichlet();
     else DslashLogFull();
-    // acceleratorCopySynchronise() is in the StencilSendToRecvFromComplete
-    //    accelerator_barrier(); 
+    acceleratorCopySynchronise();
     _grid->StencilBarrier(); 
-#ifndef ACCELERATOR_AWARE_MPI
-#warning "Using COPY VIA HOST BUFFERS IN STENCIL"
-    for(int i=0;i<Packets.size();i++){
-      if ( Packets[i].do_recv ) {
-	acceleratorCopyToDevice(Packets[i].host_recv_buf, Packets[i].recv_buf,Packets[i].rbytes);
-      }
-    }
-    _grid->HostBufferFreeAll();
-#endif
-    // run any checksums
-    for(int i=0;i<Packets.size();i++){
-      if ( Packets[i].do_recv )
-	FlightRecorder::recvLog(Packets[i].recv_buf,Packets[i].rbytes,Packets[i].from_rank);
-    }
   }
   ////////////////////////////////////////////////////////////////////////
   // Blocking send and receive. Either sequential or parallel.
@@ -502,7 +528,6 @@ public:
   template<class compressor>
   void HaloGather(const Lattice<vobj> &source,compressor &compress)
   {
-    //    accelerator_barrier();
     _grid->StencilBarrier();// Synch shared memory on a single nodes
 
     assert(source.Grid()==_grid);
@@ -515,9 +540,10 @@ public:
       compress.Point(point);
       HaloGatherDir(source,compress,point,face_idx);
     }
-    accelerator_barrier(); // All my local gathers are complete
+    accelerator_barrier();
     face_table_computed=1;
     assert(u_comm_offset==_unified_buffer_size);
+
   }
 
   /////////////////////////
@@ -553,7 +579,6 @@ public:
       accelerator_forNB(j, words, cobj::Nsimd(), {
 	  coalescedWrite(to[j] ,coalescedRead(from [j]));
       });
-      acceleratorFenceComputeStream();
     }
   }
   
@@ -644,7 +669,6 @@ public:
     for(int i=0;i<dd.size();i++){
       decompressor::DecompressFace(decompress,dd[i]);
     }
-    acceleratorFenceComputeStream(); // dependent kernels
   }
   ////////////////////////////////////////
   // Set up routines
@@ -1200,6 +1224,7 @@ public:
 	  ///////////////////////////////////////////////////////////
 	  int do_send = (comms_send|comms_partial_send) && (!shm_send );
 	  int do_recv = (comms_send|comms_partial_send) && (!shm_recv );
+	  
 	  AddPacket((void *)&send_buf[comm_off],
 		    (void *)&recv_buf[comm_off],
 		    xmit_to_rank, do_send,

Grid/tensors/Tensor_traits.h

@@ -404,12 +404,3 @@ NAMESPACE_BEGIN(Grid);
   };
 NAMESPACE_END(Grid);
 
-
-#ifdef GRID_SYCL
-template<typename T> struct
-sycl::is_device_copyable<T, typename std::enable_if<
-			      Grid::isGridTensor<T>::value  && (!std::is_trivially_copyable<T>::value),
-			      void>::type>
-  : public std::true_type {};
-#endif
-

Grid/threads/Accelerator.h

@@ -255,13 +255,17 @@ inline int  acceleratorIsCommunicable(void *ptr)
 #define GRID_SYCL_LEVEL_ZERO_IPC
 
 NAMESPACE_END(Grid);
-
-// Force deterministic reductions
-#define SYCL_REDUCTION_DETERMINISTIC
+#if 0
+#include <CL/sycl.hpp>
+#include <CL/sycl/usm.hpp>
+#include <level_zero/ze_api.h>
+#include <CL/sycl/backend/level_zero.hpp>
+#else
 #include <sycl/CL/sycl.hpp>
 #include <sycl/usm.hpp>
 #include <level_zero/ze_api.h>
 #include <sycl/ext/oneapi/backend/level_zero.hpp>
+#endif
 
 NAMESPACE_BEGIN(Grid);
 

Grid/util/FlightRecorder.cc

@@ -1,339 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid
-
-    Source file: ./lib/Init.cc
-
-    Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: Peter Boyle <peterboyle@MacBook-Pro.local>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Grid/Grid.h>
-
-NAMESPACE_BEGIN(Grid);
-///////////////////////////////////////////////////////
-// Grid Norm logging for repro testing
-///////////////////////////////////////////////////////
-int FlightRecorder::PrintEntireLog;
-int FlightRecorder::ContinueOnFail;
-int FlightRecorder::LoggingMode;
-int FlightRecorder::ChecksumComms;
-int FlightRecorder::ChecksumCommsSend;
-int32_t  FlightRecorder::XmitLoggingCounter;
-int32_t  FlightRecorder::RecvLoggingCounter;
-int32_t  FlightRecorder::CsumLoggingCounter;
-int32_t  FlightRecorder::NormLoggingCounter;
-int32_t  FlightRecorder::ReductionLoggingCounter;
-uint64_t FlightRecorder::ErrorCounter;
-std::vector<double> FlightRecorder::NormLogVector;
-std::vector<double> FlightRecorder::ReductionLogVector;
-std::vector<uint64_t> FlightRecorder::CsumLogVector;
-std::vector<uint64_t> FlightRecorder::XmitLogVector;
-std::vector<uint64_t> FlightRecorder::RecvLogVector;
-
-void FlightRecorder::ResetCounters(void)
-{
-  XmitLoggingCounter=0;
-  RecvLoggingCounter=0;
-  CsumLoggingCounter=0;
-  NormLoggingCounter=0;
-  ReductionLoggingCounter=0;
-}
-void FlightRecorder::Truncate(void)
-{
-  ResetCounters();
-  XmitLogVector.resize(0);
-  RecvLogVector.resize(0);
-  NormLogVector.resize(0);
-  CsumLogVector.resize(0);
-  ReductionLogVector.resize(0);
-}
-void FlightRecorder::SetLoggingMode(FlightRecorder::LoggingMode_t mode)
-{
-  switch ( mode ) {
-  case LoggingModePrint:
-    SetLoggingModePrint();
-    break;
-  case LoggingModeRecord:
-    SetLoggingModeRecord();
-    break;
-  case LoggingModeVerify:
-    SetLoggingModeVerify();
-    break;
-  case LoggingModeNone:
-    LoggingMode = mode;
-    Truncate();
-    break;
-  default:
-    assert(0);
-  }
-}
-
-void FlightRecorder::SetLoggingModePrint(void)
-{
-  std::cout << " FlightRecorder: set to print output " <<std::endl;
-  Truncate();
-  LoggingMode = LoggingModePrint;
-}
-void FlightRecorder::SetLoggingModeRecord(void)
-{
-  std::cout << " FlightRecorder: set to RECORD " <<std::endl;
-  Truncate();
-  LoggingMode = LoggingModeRecord;
-}
-void FlightRecorder::SetLoggingModeVerify(void)
-{
-  std::cout << " FlightRecorder: set to VERIFY " << NormLogVector.size()<< " log entries "<<std::endl;
-  ResetCounters();
-  LoggingMode = LoggingModeVerify;
-}
-uint64_t FlightRecorder::ErrorCount(void)
-{
-  return ErrorCounter;
-}
-void FlightRecorder::NormLog(double value)
-{
-  uint64_t hex = * ( (uint64_t *)&value );
-  if(LoggingMode == LoggingModePrint) {
-    std::cerr<<"FlightRecorder::NormLog : "<< NormLoggingCounter <<" "<<std::hex<< hex<<std::dec <<std::endl;
-    NormLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeRecord) {
-    std::cerr<<"FlightRecorder::NormLog RECORDING : "<< NormLoggingCounter <<" "<<std::hex<< hex<<std::dec <<std::endl;
-    NormLogVector.push_back(value);
-    NormLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeVerify) {
-
-    if(NormLoggingCounter < NormLogVector.size()){
-      uint64_t hexref  = * ( (uint64_t *)&NormLogVector[NormLoggingCounter] );
-
-      if ( (value != NormLogVector[NormLoggingCounter]) || std::isnan(value) ) {
-
-	std::cerr<<"FlightRecorder::NormLog Oops, I did it again "<< NormLoggingCounter
-		 <<std::hex<<" "<<hex<<" "<<hexref<<std::dec<<" "
-		 <<std::hexfloat<<value<<" "<< NormLogVector[NormLoggingCounter]<<std::endl;
-
-	std::cerr << " Oops got norm "<< std::hexfloat<<value<<" expect "<<NormLogVector[NormLoggingCounter] <<std::endl;
-
-	fprintf(stderr,"%s:%d Oops, I did it again! Reproduce failure for norm %d/%zu %.16e expect %.16e\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		NormLoggingCounter,NormLogVector.size(),
-		value, NormLogVector[NormLoggingCounter]); fflush(stderr);
-
-	if(!ContinueOnFail)assert(0); // Force takedown of job
-	  
-	ErrorCounter++;
-      } else {
-	if ( PrintEntireLog ) { 
-	  std::cerr<<"FlightRecorder::NormLog VALID "<< NormLoggingCounter << std::hex
-		   <<" "<<hex<<" "<<hexref
-		   <<" "<<std::hexfloat<<value<<" "<< NormLogVector[NormLoggingCounter]<<std::dec<<std::endl;
-	}
-      }
-       
-    }
-    if ( NormLogVector.size()==NormLoggingCounter ) {
-      std::cout << "FlightRecorder:: Verified entire sequence of "<<NormLoggingCounter<<" norms "<<std::endl;
-    }
-    NormLoggingCounter++;
-  }
-}
-void FlightRecorder::CsumLog(uint64_t hex)
-{
-  if(LoggingMode == LoggingModePrint) {
-    std::cerr<<"FlightRecorder::CsumLog : "<< CsumLoggingCounter <<" "<<std::hex<< hex<<std::dec <<std::endl;
-    CsumLoggingCounter++;
-  }
-
-  if(LoggingMode == LoggingModeRecord) {
-    std::cerr<<"FlightRecorder::CsumLog RECORDING : "<< NormLoggingCounter <<" "<<std::hex<< hex<<std::dec <<std::endl;
-    CsumLogVector.push_back(hex);
-    CsumLoggingCounter++;
-  }
-
-  if(LoggingMode == LoggingModeVerify) {
-    
-    if(CsumLoggingCounter < CsumLogVector.size()) {
-
-      uint64_t hexref  = CsumLogVector[CsumLoggingCounter] ;
-
-      if ( hex != hexref ) {
-
-        std::cerr<<"FlightRecorder::CsumLog Oops, I did it again "<< CsumLoggingCounter
-		 <<std::hex<<" "<<hex<<" "<<hexref<<std::dec<<std::endl;
-
-	fprintf(stderr,"%s:%d Oops, I did it again! Reproduce failure for csum %d %lx expect %lx\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		CsumLoggingCounter,hex, hexref);
-	fflush(stderr);
-
-	if(!ContinueOnFail) assert(0); // Force takedown of job
-	  
-	ErrorCounter++;
-
-      } else {
-
-	if ( PrintEntireLog ) { 
-	  std::cerr<<"FlightRecorder::CsumLog VALID "<< CsumLoggingCounter << std::hex
-		   <<" "<<hex<<" "<<hexref<<std::dec<<std::endl;
-	}
-      }
-    }  
-    if ( CsumLogVector.size()==CsumLoggingCounter ) {
-      std::cout << "FlightRecorder:: Verified entire sequence of "<<CsumLoggingCounter<<" checksums "<<std::endl;
-    }
-    CsumLoggingCounter++;
-  }
-}
-void FlightRecorder::ReductionLog(double local,double global)
-{
-  uint64_t hex_l = * ( (uint64_t *)&local );
-  uint64_t hex_g = * ( (uint64_t *)&global );
-  if(LoggingMode == LoggingModePrint) {
-    std::cerr<<"FlightRecorder::ReductionLog : "<< ReductionLoggingCounter <<" "<< std::hex << hex_l << " -> " <<hex_g<<std::dec <<std::endl;
-    ReductionLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeRecord) {
-    std::cerr<<"FlightRecorder::ReductionLog RECORDING : "<< ReductionLoggingCounter <<" "<< std::hex << hex_l << " -> " <<hex_g<<std::dec <<std::endl;
-    ReductionLogVector.push_back(global);
-    ReductionLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeVerify) {
-    if(ReductionLoggingCounter < ReductionLogVector.size()){
-      if ( global != ReductionLogVector[ReductionLoggingCounter] ) {
-	fprintf(stderr,"%s:%d Oops, MPI_Allreduce did it again! Reproduce failure for norm %d/%zu glb %.16e lcl %.16e expect glb %.16e\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		ReductionLoggingCounter,ReductionLogVector.size(),
-		global, local, ReductionLogVector[ReductionLoggingCounter]); fflush(stderr);
-	
-	if ( !ContinueOnFail ) assert(0);
-
-	ErrorCounter++;
-      } else {
-	if ( PrintEntireLog ) { 
-	  std::cerr<<"FlightRecorder::ReductionLog : VALID "<< ReductionLoggingCounter <<" "<< std::hexfloat << local << "-> "<< global <<std::endl;
-	}
-      }
-    }
-    if ( ReductionLogVector.size()==ReductionLoggingCounter ) {
-      std::cout << "FlightRecorder::ReductionLog : Verified entire sequence of "<<ReductionLoggingCounter<<" norms "<<std::endl;
-    }
-    ReductionLoggingCounter++;
-  }
-}
-void FlightRecorder::xmitLog(void *buf,uint64_t bytes)
-{
-  if ( ChecksumCommsSend ){
-  uint64_t *ubuf = (uint64_t *)buf;
-  if(LoggingMode == LoggingModeNone) return;
-#ifdef GRID_SYCL
-  uint64_t _xor = svm_xor(ubuf,bytes/sizeof(uint64_t));
-  if(LoggingMode == LoggingModePrint) {
-    std::cerr<<"FlightRecorder::xmitLog : "<< XmitLoggingCounter <<" "<< std::hex << _xor <<std::dec <<std::endl;
-    XmitLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeRecord) {
-    std::cerr<<"FlightRecorder::xmitLog RECORD : "<< XmitLoggingCounter <<" "<< std::hex << _xor <<std::dec <<std::endl;
-    XmitLogVector.push_back(_xor);
-    XmitLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeVerify) {
-    if(XmitLoggingCounter < XmitLogVector.size()){
-      if ( _xor != XmitLogVector[XmitLoggingCounter] ) {
-	fprintf(stderr,"%s:%d Oops, send buf difference! Reproduce failure for xmit %d/%zu  %lx expect glb %lx\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		XmitLoggingCounter,XmitLogVector.size(),
-		_xor, XmitLogVector[XmitLoggingCounter]); fflush(stderr);
-	
-	if ( !ContinueOnFail ) assert(0);
-
-	ErrorCounter++;
-      } else {
-	if ( PrintEntireLog ) { 
-	  std::cerr<<"FlightRecorder::XmitLog : VALID "<< XmitLoggingCounter <<" "<< std::hexfloat << _xor << " "<<  XmitLogVector[XmitLoggingCounter] <<std::endl;
-	}
-      }
-    }
-    if ( XmitLogVector.size()==XmitLoggingCounter ) {
-      std::cout << "FlightRecorder::ReductionLog : Verified entire sequence of "<<XmitLoggingCounter<<" sends "<<std::endl;
-    }
-    XmitLoggingCounter++;
-  }
-#endif
-  } else {
-    uint64_t word = 1;
-    deviceVector<uint64_t> dev(1);
-    acceleratorCopyToDevice(&word,&dev[0],sizeof(uint64_t));
-    acceleratorCopySynchronise();
-    MPI_Barrier(MPI_COMM_WORLD);
-  }
-}
-void FlightRecorder::recvLog(void *buf,uint64_t bytes,int rank)
-{
-  if ( ChecksumComms ){
-  uint64_t *ubuf = (uint64_t *)buf;
-  if(LoggingMode == LoggingModeNone) return;
-#ifdef GRID_SYCL
-  uint64_t _xor = svm_xor(ubuf,bytes/sizeof(uint64_t));
-  if(LoggingMode == LoggingModePrint) {
-    std::cerr<<"FlightRecorder::recvLog : "<< RecvLoggingCounter <<" "<< std::hex << _xor <<std::dec <<std::endl;
-    RecvLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeRecord) {
-    std::cerr<<"FlightRecorder::recvLog RECORD : "<< RecvLoggingCounter <<" "<< std::hex << _xor <<std::dec <<std::endl;
-    RecvLogVector.push_back(_xor);
-    RecvLoggingCounter++;
-  }
-  if(LoggingMode == LoggingModeVerify) {
-    if(RecvLoggingCounter < RecvLogVector.size()){
-      if ( _xor != RecvLogVector[RecvLoggingCounter] ) {
-	fprintf(stderr,"%s:%d Oops, recv buf difference! Reproduce failure for recv %d/%zu  %lx expect glb %lx from MPI rank %d\n",
-		GridHostname(),
-		GlobalSharedMemory::WorldShmRank,
-		RecvLoggingCounter,RecvLogVector.size(),
-		_xor, RecvLogVector[RecvLoggingCounter],rank); fflush(stderr);
-	
-	if ( !ContinueOnFail ) assert(0);
-
-	ErrorCounter++;
-      } else {
-	if ( PrintEntireLog ) { 
-	  std::cerr<<"FlightRecorder::RecvLog : VALID "<< RecvLoggingCounter <<" "<< std::hexfloat << _xor << " "<<  RecvLogVector[RecvLoggingCounter] <<std::endl;
-	}
-      }
-    }
-    if ( RecvLogVector.size()==RecvLoggingCounter ) {
-      std::cout << "FlightRecorder::ReductionLog : Verified entire sequence of "<<RecvLoggingCounter<<" sends "<<std::endl;
-    }
-    RecvLoggingCounter++;
-  }
-#endif
-  }
-}
-
-NAMESPACE_END(Grid);

Grid/util/FlightRecorder.h

@@ -1,43 +0,0 @@
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-class FlightRecorder {
- public:
-  enum LoggingMode_t {
-    LoggingModeNone,
-    LoggingModePrint,
-    LoggingModeRecord,
-    LoggingModeVerify
-  };
-  
-  static int                   LoggingMode;
-  static uint64_t              ErrorCounter;
-  static int32_t               XmitLoggingCounter;
-  static int32_t               RecvLoggingCounter;
-  static int32_t               CsumLoggingCounter;
-  static int32_t               NormLoggingCounter;
-  static int32_t               ReductionLoggingCounter;
-  static std::vector<uint64_t> XmitLogVector;
-  static std::vector<uint64_t> RecvLogVector;
-  static std::vector<uint64_t> CsumLogVector;
-  static std::vector<double>   NormLogVector;
-  static std::vector<double>   ReductionLogVector;
-  static int ContinueOnFail;
-  static int PrintEntireLog;
-  static int ChecksumComms;
-  static int ChecksumCommsSend;
-  static void SetLoggingModePrint(void);
-  static void SetLoggingModeRecord(void);
-  static void SetLoggingModeVerify(void);
-  static void SetLoggingMode(LoggingMode_t mode);
-  static void NormLog(double value);
-  static void CsumLog(uint64_t csum);
-  static void ReductionLog(double lcl, double glbl);
-  static void Truncate(void);
-  static void ResetCounters(void);
-  static uint64_t ErrorCount(void);
-  static void xmitLog(void *,uint64_t bytes);
-  static void recvLog(void *,uint64_t bytes,int rank);
-};
-NAMESPACE_END(Grid);
-

Grid/util/Init.cc

@@ -77,10 +77,6 @@ feenableexcept (unsigned int excepts)
 }
 #endif
 
-#ifndef HOST_NAME_MAX
-#define HOST_NAME_MAX _POSIX_HOST_NAME_MAX
-#endif
-
 NAMESPACE_BEGIN(Grid);
 
 //////////////////////////////////////////////////////
@@ -94,12 +90,7 @@ int GridThread::_threads =1;
 int GridThread::_hyperthreads=1;
 int GridThread::_cores=1;
 
-char hostname[HOST_NAME_MAX+1];
 
-char *GridHostname(void)
-{
-  return hostname;
-}
 const Coordinate &GridDefaultLatt(void)     {return Grid_default_latt;};
 const Coordinate &GridDefaultMpi(void)      {return Grid_default_mpi;};
 const Coordinate GridDefaultSimd(int dims,int nsimd)
@@ -402,8 +393,6 @@ void Grid_init(int *argc,char ***argv)
   std::cout << GridLogMessage << "MPI is initialised and logging filters activated "<<std::endl;
   std::cout << GridLogMessage << "================================================ "<<std::endl;
 
-  gethostname(hostname, HOST_NAME_MAX+1);
-  std::cout << GridLogMessage << "This rank is running on host "<< hostname<<std::endl;
 
   /////////////////////////////////////////////////////////
   // Reporting

Grid/util/Init.h

@@ -34,8 +34,6 @@ NAMESPACE_BEGIN(Grid);
 void Grid_init(int *argc,char ***argv);
 void Grid_finalize(void);
 
-char * GridHostname(void);
-
 // internal, controled with --handle
 void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr);
 void Grid_debug_handler_init(void);
@@ -70,6 +68,5 @@ void GridParseLayout(char **argv,int argc,
 void printHash(void);
 
 
-
 NAMESPACE_END(Grid);
 

Grid/util/Util.h

@@ -1,6 +1,6 @@
-#pragma once
+#ifndef GRID_UTIL_H
+#define GRID_UTIL_H
 #include <Grid/util/Coordinate.h>
 #include <Grid/util/Lexicographic.h>
 #include <Grid/util/Init.h>
-#include <Grid/util/FlightRecorder.h>
-
+#endif

benchmarks/Benchmark_usqcd.cc

@@ -1,968 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./benchmarks/Benchmark_usqcd.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: paboyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#include <Grid/Grid.h>
-#include <Grid/algorithms/blas/BatchedBlas.h>
-
-using namespace Grid;
-
-std::vector<int> L_list;
-std::vector<int> Ls_list;
-std::vector<double> mflop_list;
-
-double mflop_ref;
-double mflop_ref_err;
-
-int NN_global;
-
-FILE * FP;
-
-struct time_statistics{
-  double mean;
-  double err;
-  double min;
-  double max;
-
-  void statistics(std::vector<double> v){
-      double sum = std::accumulate(v.begin(), v.end(), 0.0);
-      mean = sum / v.size();
-
-      std::vector<double> diff(v.size());
-      std::transform(v.begin(), v.end(), diff.begin(), [=](double x) { return x - mean; });
-      double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
-      err = std::sqrt(sq_sum / (v.size()*(v.size() - 1)));
-
-      auto result = std::minmax_element(v.begin(), v.end());
-      min = *result.first;
-      max = *result.second;
-}
-};
-
-void comms_header(){
-  std::cout <<GridLogMessage << " L  "<<"\t"<<" Ls  "<<"\t"
-            <<"bytes\t MB/s uni  \t\t MB/s bidi "<<std::endl;
-};
-
-struct controls {
-  int Opt;
-  int CommsOverlap;
-  Grid::CartesianCommunicator::CommunicatorPolicy_t CommsAsynch;
-};
-
-class Benchmark {
-public:
-  static void Decomposition (void ) {
-
-    int threads = GridThread::GetThreads();
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "= Grid is setup to use "<<threads<<" threads"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage<<"Grid Default Decomposition patterns\n";
-    std::cout<<GridLogMessage<<"\tOpenMP threads : "<<GridThread::GetThreads()<<std::endl;
-    std::cout<<GridLogMessage<<"\tMPI tasks      : "<<GridCmdVectorIntToString(GridDefaultMpi())<<std::endl;
-    std::cout<<GridLogMessage<<"\tvReal          : "<<sizeof(vReal )*8    <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vReal::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage<<"\tvRealF         : "<<sizeof(vRealF)*8    <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealF::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage<<"\tvRealD         : "<<sizeof(vRealD)*8    <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealD::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage<<"\tvComplex       : "<<sizeof(vComplex )*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplex::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage<<"\tvComplexF      : "<<sizeof(vComplexF)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexF::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage<<"\tvComplexD      : "<<sizeof(vComplexD)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexD::Nsimd()))<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-  }
-
-  static void Comms(void)
-  {
-    int Nloop=200;
-    int nmu=0;
-    int maxlat=32;
-
-    Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
-    Coordinate mpi_layout  = GridDefaultMpi();
-
-    for(int mu=0;mu<Nd;mu++) if (mpi_layout[mu]>1) nmu++;
-
-    std::vector<double> t_time(Nloop);
-    time_statistics timestat;
-
-    std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "= Benchmarking threaded STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
-    std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
-    comms_header();
-
-    fprintf(FP,"Communications\n\n");
-    fprintf(FP,"Packet bytes, direction, GB/s per node\n");
-    for(int lat=16;lat<=maxlat;lat+=8){
-      //      for(int Ls=8;Ls<=8;Ls*=2){
-      { int Ls=12;
-
-	Coordinate latt_size  ({lat*mpi_layout[0],
-	      lat*mpi_layout[1],
-	      lat*mpi_layout[2],
-	      lat*mpi_layout[3]});
-
-	GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-	RealD Nrank = Grid._Nprocessors;
-	RealD Nnode = Grid.NodeCount();
-	RealD ppn = Nrank/Nnode;
-
-	std::vector<HalfSpinColourVectorD *> xbuf(8);
-	std::vector<HalfSpinColourVectorD *> rbuf(8);
-	//Grid.ShmBufferFreeAll();
-	uint64_t bytes=lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD);
-	for(int d=0;d<8;d++){
-	  xbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes);
-	  rbuf[d] = (HalfSpinColourVectorD *)acceleratorAllocDevice(bytes);
-	  //	  bzero((void *)xbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	  //	  bzero((void *)rbuf[d],lat*lat*lat*Ls*sizeof(HalfSpinColourVectorD));
-	}
-
-	//	int ncomm;
-	double dbytes;
-
-        for(int dir=0;dir<8;dir++) {
-	  int mu =dir % 4;
-	  if (mpi_layout[mu]>1 ) {
-
-	    std::vector<double> times(Nloop);
-	    for(int i=0;i<Nloop;i++){
-
-	      dbytes=0;	        
-	      double start=usecond();
-	      int xmit_to_rank;
-	      int recv_from_rank;
-
-	      if ( dir == mu ) { 
-		int comm_proc=1;
-		Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	      } else { 
-		int comm_proc = mpi_layout[mu]-1;
-		Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
-	      }
-	      Grid.SendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank,
-				  (void *)&rbuf[dir][0], recv_from_rank,
-				  bytes);
-	      dbytes+=bytes;
-	     
-	      double stop=usecond();
-	      t_time[i] = stop-start; // microseconds
-
-	    }
-	    timestat.statistics(t_time);
-	  
-	    dbytes=dbytes*ppn;
-	    double xbytes    = dbytes*0.5;
-	    double bidibytes = dbytes;
-	  
-	    std::cout<<GridLogMessage << lat<<"\t"<<Ls<<"\t "
-		     << bytes << " \t "
-		     <<xbytes/timestat.mean
-		     << "\t\t"
-		     << bidibytes/timestat.mean<< std::endl;
-	    fprintf(FP,"%ld, %d, %f\n",(long)bytes,dir,bidibytes/timestat.mean/1000.);
-	  }
-	}
-	for(int d=0;d<8;d++){
-	  acceleratorFreeDevice(xbuf[d]);
-	  acceleratorFreeDevice(rbuf[d]);
-	}
-      }
-    }
-    fprintf(FP,"\n\n");
-    
-    return;
-  }
-
-  
-  static void Memory(void)
-  {
-    const int Nvec=8;
-    typedef Lattice< iVector< vReal,Nvec> > LatticeVec;
-    typedef iVector<vReal,Nvec> Vec;
-
-    Coordinate simd_layout = GridDefaultSimd(Nd,vReal::Nsimd());
-    Coordinate mpi_layout  = GridDefaultMpi();
-
-    fprintf(FP,"Memory Bandwidth\n\n");
-    fprintf(FP,"Bytes, GB/s per node\n");
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "= Benchmarking a*x + y bandwidth"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<"\t\t seconds"<< "\t\tGB/s / node"<<std::endl;
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-  
-    //    uint64_t NP;
-    uint64_t NN;
-
-
-  uint64_t lmax=40;
-#define NLOOP (1000*lmax*lmax*lmax*lmax/lat/lat/lat/lat)
-
-    GridSerialRNG          sRNG;      sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-    for(int lat=8;lat<=lmax;lat+=8){
-
-      Coordinate latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
-      int64_t vol= latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
-
-      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-
-      //      NP= Grid.RankCount();
-      NN =Grid.NodeCount();
-
-      Vec rn ; random(sRNG,rn);
-
-      LatticeVec z(&Grid); z=Zero();
-      LatticeVec x(&Grid); x=Zero();
-      LatticeVec y(&Grid); y=Zero();
-      double a=2.0;
-
-      uint64_t Nloop=NLOOP;
-
-      double start=usecond();
-      for(int i=0;i<Nloop;i++){
-	z=a*x-y;
-      }
-      double stop=usecond();
-      double time = (stop-start)/Nloop*1000;
-     
-      double flops=vol*Nvec*2;// mul,add
-      double bytes=3.0*vol*Nvec*sizeof(Real);
-      std::cout<<GridLogMessage<<std::setprecision(3) 
-	       << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t"<<flops/time<<"\t\t"<<(stop-start)/1000./1000.
-	       << "\t\t"<< bytes/time/NN <<std::endl;
-
-      fprintf(FP,"%ld, %f\n",(long)bytes,bytes/time/NN);
-
-    }
-    fprintf(FP,"\n\n");
-  };
-
-
-  static void BLAS(void)
-  {
-    //int nbasis, int nrhs, int coarseVol
-    int  basis[] = { 16,32,64 };
-    int  rhs[]   = { 8,16,32 };
-    int  vol  = 4*4*4*4;
-
-    GridBLAS blas;
-    
-    std::cout<<GridLogMessage << "=================================================================================="<<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\n");
-
-    for(int b=0;b<3;b++){
-    for(int r=0;r<3;r++){
-      int M=basis[b];
-      int N=rhs[r];
-      int K=basis[b];
-      int BATCH=vol;
-      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) 
-	       << M<<"\t\t"<<N<<"\t\t"<<K<<"\t\t"<<BATCH<<"\t\t"<<p<<std::endl;
-    }}
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-    std::cout<<GridLogMessage << "  M  "<<"\t\t"<<"N"<<"\t\t\t"<<"K"<<"\t\t"<<"Gflop/s / rank (block project)"<<std::endl;
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-    for(int b=0;b<3;b++){
-    for(int r=0;r<3;r++){
-      int M=basis[b];
-      int N=rhs[r];
-      int K=vol;
-      int BATCH=vol;
-      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) 
-	       << M<<"\t\t"<<N<<"\t\t"<<K<<"\t\t"<<BATCH<<"\t\t"<<p<<std::endl;
-    }}
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-    std::cout<<GridLogMessage << "  M  "<<"\t\t"<<"N"<<"\t\t\t"<<"K"<<"\t\t"<<"Gflop/s / rank (block promote)"<<std::endl;
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-    for(int b=0;b<3;b++){
-    for(int r=0;r<3;r++){
-      int M=rhs[r];
-      int N=vol;
-      int K=basis[b];
-      int BATCH=vol;
-      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) 
-	       << M<<"\t\t"<<N<<"\t\t"<<K<<"\t\t"<<BATCH<<"\t\t"<<p<<std::endl;
-    }}
-    fprintf(FP,"\n\n\n");
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  };
-  
-
-  static void SU4(void)
-  {
-    const int Nc4=4;
-    typedef Lattice< iMatrix< vComplexF,Nc4> > LatticeSU4;
-
-    Coordinate simd_layout = GridDefaultSimd(Nd,vComplexF::Nsimd());
-    Coordinate mpi_layout  = GridDefaultMpi();
-    
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "= Benchmarking z = y*x SU(4) bandwidth"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<"\t\t seconds"<< "\t\tGB/s / node"<<std::endl;
-    std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-  
-    uint64_t NN;
-
-
-    uint64_t lmax=32;
-
-    GridSerialRNG          sRNG;      sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-    for(int lat=8;lat<=lmax;lat+=8){
-
-      Coordinate latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
-      int64_t vol= latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
-
-      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
-
-      NN =Grid.NodeCount();
-
-
-      LatticeSU4 z(&Grid); z=Zero();
-      LatticeSU4 x(&Grid); x=Zero();
-      LatticeSU4 y(&Grid); y=Zero();
-      //      double a=2.0;
-
-      uint64_t Nloop=NLOOP;
-
-      double start=usecond();
-      for(int i=0;i<Nloop;i++){
-	z=x*y;
-      }
-      double stop=usecond();
-      double time = (stop-start)/Nloop*1000;
-     
-      double flops=vol*Nc4*Nc4*(6+(Nc4-1)*8);// mul,add
-      double bytes=3.0*vol*Nc4*Nc4*2*sizeof(RealF);
-      std::cout<<GridLogMessage<<std::setprecision(3) 
-	       << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t"<<flops/time<<"\t\t"<<(stop-start)/1000./1000.
-	       << "\t\t"<< bytes/time/NN <<std::endl;
-
-    }
-  };
-
-
-  static double DWF(int Ls,int L)
-  {
-    RealD mass=0.1;
-    RealD M5  =1.8;
-
-    double mflops;
-    double mflops_best = 0;
-    double mflops_worst= 0;
-    std::vector<double> mflops_all;
-
-    ///////////////////////////////////////////////////////
-    // Set/Get the layout & grid size
-    ///////////////////////////////////////////////////////
-    int threads = GridThread::GetThreads();
-    Coordinate mpi = GridDefaultMpi(); assert(mpi.size()==4);
-    Coordinate local({L,L,L,L});
-    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
-
-    GridCartesian         * TmpGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, 
-								       GridDefaultSimd(Nd,vComplex::Nsimd()),
-								       GridDefaultMpi());
-    uint64_t NP = TmpGrid->RankCount();
-    uint64_t NN = TmpGrid->NodeCount();
-    NN_global=NN;
-    uint64_t SHM=NP/NN;
-
-
-    ///////// Welcome message ////////////
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "Benchmark DWF on "<<L<<"^4 local volume "<<std::endl;
-    std::cout<<GridLogMessage << "* Nc             : "<<Nc<<std::endl;
-    std::cout<<GridLogMessage << "* Global volume  : "<<GridCmdVectorIntToString(latt4)<<std::endl;
-    std::cout<<GridLogMessage << "* Ls             : "<<Ls<<std::endl;
-    std::cout<<GridLogMessage << "* ranks          : "<<NP  <<std::endl;
-    std::cout<<GridLogMessage << "* nodes          : "<<NN  <<std::endl;
-    std::cout<<GridLogMessage << "* ranks/node     : "<<SHM <<std::endl;
-    std::cout<<GridLogMessage << "* ranks geom     : "<<GridCmdVectorIntToString(mpi)<<std::endl;
-    std::cout<<GridLogMessage << "* Using "<<threads<<" threads"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    ///////// Lattice Init ////////////
-    GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-    GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-    
-    ///////// RNG Init ////////////
-    std::vector<int> seeds4({1,2,3,4});
-    std::vector<int> seeds5({5,6,7,8});
-    GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
-    GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
-    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
-
-    typedef DomainWallFermionF Action;
-    typedef typename Action::FermionField Fermion;
-    typedef LatticeGaugeFieldF Gauge;
-    
-    ///////// Source preparation ////////////
-    Gauge Umu(UGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
-    Fermion src   (FGrid); random(RNG5,src);
-    Fermion src_e (FrbGrid);
-    Fermion src_o (FrbGrid);
-    Fermion r_e   (FrbGrid);
-    Fermion r_o   (FrbGrid);
-    Fermion r_eo  (FGrid);
-    Action Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-
-    {
-
-      pickCheckerboard(Even,src_e,src);
-      pickCheckerboard(Odd,src_o,src);
-
-#ifdef AVX512
-      const int num_cases = 3;
-#else 
-      const int num_cases = 2;
-#endif      
-      std::string fmt("G/S/C ; G/O/C ; G/S/S ; G/O/S ");
-
-      controls Cases [] = {
-	{  WilsonKernelsStatic::OptGeneric   ,  WilsonKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent },
-	{  WilsonKernelsStatic::OptHandUnroll,  WilsonKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent },
-	{  WilsonKernelsStatic::OptInlineAsm ,  WilsonKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent }
-      }; 
-
-      for(int c=0;c<num_cases;c++) {
-	
-	WilsonKernelsStatic::Comms = Cases[c].CommsOverlap;
-	WilsonKernelsStatic::Opt   = Cases[c].Opt;
-	CartesianCommunicator::SetCommunicatorPolicy(Cases[c].CommsAsynch);
-
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
-	if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using ASM      WilsonKernels" <<std::endl;
-	if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) std::cout << GridLogMessage<< "* Using UNROLLED WilsonKernels" <<std::endl;
-	if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
-	if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential Comms/Compute" <<std::endl;
-	std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-	int nwarm = 10;
-	double t0=usecond();
-	FGrid->Barrier();
-	for(int i=0;i<nwarm;i++){
-	  Dw.DhopEO(src_o,r_e,DaggerNo);
-	}
-	FGrid->Barrier();
-	double t1=usecond();
-	uint64_t ncall = 500;
-
-	FGrid->Broadcast(0,&ncall,sizeof(ncall));
-
-	//	std::cout << GridLogMessage << " Estimate " << ncall << " calls per second"<<std::endl;
-
-	time_statistics timestat;
-	std::vector<double> t_time(ncall);
-	for(uint64_t i=0;i<ncall;i++){
-	  t0=usecond();
-	  Dw.DhopEO(src_o,r_e,DaggerNo);
-	  t1=usecond();
-	  t_time[i] = t1-t0;
-	}
-	FGrid->Barrier();
-	
-	double volume=Ls;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
-
-	// Nc=3 gives
-	// 1344= 3*(2*8+6)*2*8 + 8*3*2*2 + 3*4*2*8
-	// 1344 = Nc* (6+(Nc-1)*8)*2*Nd + Nd*Nc*2*2  + Nd*Nc*Ns*2
-	//	double flops=(1344.0*volume)/2;
-	double fps = Nc* (6+(Nc-1)*8)*Ns*Nd + 2*Nd*Nc*Ns  + 2*Nd*Nc*Ns*2;
-
-	double flops=(fps*volume)/2;
-	double mf_hi, mf_lo, mf_err;
-
-	timestat.statistics(t_time);
-	mf_hi = flops/timestat.min;
-	mf_lo = flops/timestat.max;
-	mf_err= flops/timestat.min * timestat.err/timestat.mean;
-
-	mflops = flops/timestat.mean;
-	mflops_all.push_back(mflops);
-	if ( mflops_best == 0   ) mflops_best = mflops;
-	if ( mflops_worst== 0   ) mflops_worst= mflops;
-	if ( mflops>mflops_best ) mflops_best = mflops;
-	if ( mflops<mflops_worst) mflops_worst= mflops;
-
-	std::cout<<GridLogMessage<< "Deo FlopsPerSite is "<<fps<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s =   "<< mflops << " ("<<mf_err<<") " << mf_lo<<"-"<<mf_hi <<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per rank   "<< mflops/NP<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per node   "<< mflops/NN<<std::endl;
-
-      }
-
-      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-      std::cout<<GridLogMessage << L<<"^4 x "<<Ls<< " Deo Best  mflop/s        =   "<< mflops_best << " ; " << mflops_best/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage << L<<"^4 x "<<Ls<< " Deo Worst mflop/s        =   "<< mflops_worst<< " ; " << mflops_worst/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage <<fmt << std::endl;
-      std::cout<<GridLogMessage ;
-
-      for(int i=0;i<mflops_all.size();i++){
-	std::cout<<mflops_all[i]/NN<<" ; " ;
-      }
-      std::cout<<std::endl;
-      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    }
-    return mflops_best;
-  }
-
-
-  static double Staggered(int L)
-  {
-    double mflops;
-    double mflops_best = 0;
-    double mflops_worst= 0;
-    std::vector<double> mflops_all;
-
-    ///////////////////////////////////////////////////////
-    // Set/Get the layout & grid size
-    ///////////////////////////////////////////////////////
-    int threads = GridThread::GetThreads();
-    Coordinate mpi = GridDefaultMpi(); assert(mpi.size()==4);
-    Coordinate local({L,L,L,L});
-    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
-    
-    GridCartesian         * TmpGrid   = SpaceTimeGrid::makeFourDimGrid(latt4,
-								       GridDefaultSimd(Nd,vComplex::Nsimd()),
-								       GridDefaultMpi());
-    uint64_t NP = TmpGrid->RankCount();
-    uint64_t NN = TmpGrid->NodeCount();
-    NN_global=NN;
-    uint64_t SHM=NP/NN;
-
-
-    ///////// Welcome message ////////////
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "Benchmark ImprovedStaggered on "<<L<<"^4 local volume "<<std::endl;
-    std::cout<<GridLogMessage << "* Global volume  : "<<GridCmdVectorIntToString(latt4)<<std::endl;
-    std::cout<<GridLogMessage << "* ranks          : "<<NP  <<std::endl;
-    std::cout<<GridLogMessage << "* nodes          : "<<NN  <<std::endl;
-    std::cout<<GridLogMessage << "* ranks/node     : "<<SHM <<std::endl;
-    std::cout<<GridLogMessage << "* ranks geom     : "<<GridCmdVectorIntToString(mpi)<<std::endl;
-    std::cout<<GridLogMessage << "* Using "<<threads<<" threads"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    ///////// Lattice Init ////////////
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
-    
-    ///////// RNG Init ////////////
-    std::vector<int> seeds4({1,2,3,4});
-    GridParallelRNG          RNG4(FGrid);  RNG4.SeedFixedIntegers(seeds4);
-    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
-
-    RealD mass=0.1;
-    RealD c1=9.0/8.0;
-    RealD c2=-1.0/24.0;
-    RealD u0=1.0;
-
-    typedef ImprovedStaggeredFermionF Action;
-    typedef typename Action::FermionField Fermion; 
-    typedef LatticeGaugeFieldF Gauge;
-    
-    Gauge Umu(FGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
-
-    typename Action::ImplParams params;
-    Action Ds(Umu,Umu,*FGrid,*FrbGrid,mass,c1,c2,u0,params);
-
-    ///////// Source preparation ////////////
-    Fermion src   (FGrid); random(RNG4,src);
-    Fermion src_e (FrbGrid);
-    Fermion src_o (FrbGrid);
-    Fermion r_e   (FrbGrid);
-    Fermion r_o   (FrbGrid);
-    Fermion r_eo  (FGrid);
-  
-    {
-
-      pickCheckerboard(Even,src_e,src);
-      pickCheckerboard(Odd,src_o,src);
-    
-      const int num_cases = 2;
-      std::string fmt("G/S/C ; G/O/C ; G/S/S ; G/O/S ");
-      
-      controls Cases [] = {
-	{  StaggeredKernelsStatic::OptGeneric   ,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
-	{  StaggeredKernelsStatic::OptHandUnroll,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
-	{  StaggeredKernelsStatic::OptInlineAsm ,  StaggeredKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  }
-      }; 
-
-      for(int c=0;c<num_cases;c++) {
-	
-	StaggeredKernelsStatic::Comms = Cases[c].CommsOverlap;
-	StaggeredKernelsStatic::Opt   = Cases[c].Opt;
-	CartesianCommunicator::SetCommunicatorPolicy(Cases[c].CommsAsynch);
-      
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	if ( StaggeredKernelsStatic::Opt == StaggeredKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc StaggeredKernels" <<std::endl;
-	std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	
-	int nwarm = 10;
-	double t0=usecond();
-	FGrid->Barrier();
-	for(int i=0;i<nwarm;i++){
-	  Ds.DhopEO(src_o,r_e,DaggerNo);
-	}
-	FGrid->Barrier();
-	double t1=usecond();
-	uint64_t ncall = 500;
-
-	FGrid->Broadcast(0,&ncall,sizeof(ncall));
-
-	//	std::cout << GridLogMessage << " Estimate " << ncall << " calls per second"<<std::endl;
-
-	time_statistics timestat;
-	std::vector<double> t_time(ncall);
-	for(uint64_t i=0;i<ncall;i++){
-	  t0=usecond();
-	  Ds.DhopEO(src_o,r_e,DaggerNo);
-	  t1=usecond();
-	  t_time[i] = t1-t0;
-	}
-	FGrid->Barrier();
-	
-	double volume=1;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
-	double flops=(1146.0*volume)/2;
-	double mf_hi, mf_lo, mf_err;
-	
-	timestat.statistics(t_time);
-	mf_hi = flops/timestat.min;
-	mf_lo = flops/timestat.max;
-	mf_err= flops/timestat.min * timestat.err/timestat.mean;
-
-	mflops = flops/timestat.mean;
-	mflops_all.push_back(mflops);
-	if ( mflops_best == 0   ) mflops_best = mflops;
-	if ( mflops_worst== 0   ) mflops_worst= mflops;
-	if ( mflops>mflops_best ) mflops_best = mflops;
-	if ( mflops<mflops_worst) mflops_worst= mflops;
-	
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s =   "<< mflops << " ("<<mf_err<<") " << mf_lo<<"-"<<mf_hi <<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per rank   "<< mflops/NP<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Deo mflop/s per node   "<< mflops/NN<<std::endl;
-      
-      }
-
-      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Best  mflop/s        =   "<< mflops_best << " ; " << mflops_best/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Worst mflop/s        =   "<< mflops_worst<< " ; " << mflops_worst/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage <<fmt << std::endl;
-      std::cout<<GridLogMessage ;
-
-      for(int i=0;i<mflops_all.size();i++){
-	std::cout<<mflops_all[i]/NN<<" ; " ;
-      }
-      std::cout<<std::endl;
-    }
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    return mflops_best;
-  }
-
-  static double Clover(int L)
-  {
-    double mflops;
-    double mflops_best = 0;
-    double mflops_worst= 0;
-    std::vector<double> mflops_all;
-
-    ///////////////////////////////////////////////////////
-    // Set/Get the layout & grid size
-    ///////////////////////////////////////////////////////
-    int threads = GridThread::GetThreads();
-    Coordinate mpi = GridDefaultMpi(); assert(mpi.size()==4);
-    Coordinate local({L,L,L,L});
-    Coordinate latt4({local[0]*mpi[0],local[1]*mpi[1],local[2]*mpi[2],local[3]*mpi[3]});
-    
-    GridCartesian         * TmpGrid   = SpaceTimeGrid::makeFourDimGrid(latt4,
-								       GridDefaultSimd(Nd,vComplex::Nsimd()),
-								       GridDefaultMpi());
-    uint64_t NP = TmpGrid->RankCount();
-    uint64_t NN = TmpGrid->NodeCount();
-    NN_global=NN;
-    uint64_t SHM=NP/NN;
-
-
-    ///////// Welcome message ////////////
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << "Benchmark Clover on "<<L<<"^4 local volume "<<std::endl;
-    std::cout<<GridLogMessage << "* Global volume  : "<<GridCmdVectorIntToString(latt4)<<std::endl;
-    std::cout<<GridLogMessage << "* ranks          : "<<NP  <<std::endl;
-    std::cout<<GridLogMessage << "* nodes          : "<<NN  <<std::endl;
-    std::cout<<GridLogMessage << "* ranks/node     : "<<SHM <<std::endl;
-    std::cout<<GridLogMessage << "* ranks geom     : "<<GridCmdVectorIntToString(mpi)<<std::endl;
-    std::cout<<GridLogMessage << "* Using "<<threads<<" threads"<<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-
-    ///////// Lattice Init ////////////
-    GridCartesian         * FGrid   = SpaceTimeGrid::makeFourDimGrid(latt4, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-    GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
-    
-    ///////// RNG Init ////////////
-    std::vector<int> seeds4({1,2,3,4});
-    GridParallelRNG          RNG4(FGrid);  RNG4.SeedFixedIntegers(seeds4);
-    std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
-
-    RealD mass=0.1;
-    RealD csw=1.0;
-
-    typedef WilsonCloverFermionF Action;
-    typedef typename Action::FermionField Fermion; 
-    typedef LatticeGaugeFieldF Gauge;
-    
-    Gauge Umu(FGrid);  SU<Nc>::HotConfiguration(RNG4,Umu); 
-
-    Action Dc(Umu,*FGrid,*FrbGrid,mass,csw,csw);
-
-    ///////// Source preparation ////////////
-    Fermion src   (FGrid); random(RNG4,src);
-    Fermion r     (FGrid);
-  
-    {
-
-      const int num_cases = 1;
-      std::string fmt("G/S/C ; G/O/C ; G/S/S ; G/O/S ");
-      
-      controls Cases [] = {
-	{  WilsonKernelsStatic::OptGeneric   ,  WilsonKernelsStatic::CommsAndCompute  ,CartesianCommunicator::CommunicatorPolicyConcurrent  },
-      }; 
-
-      for(int c=0;c<num_cases;c++) {
-	
-	WilsonKernelsStatic::Comms = Cases[c].CommsOverlap;
-	WilsonKernelsStatic::Opt   = Cases[c].Opt;
-	CartesianCommunicator::SetCommunicatorPolicy(Cases[c].CommsAsynch);
-      
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
-	std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-	
-	int nwarm = 10;
-	double t0=usecond();
-	FGrid->Barrier();
-	for(int i=0;i<nwarm;i++){
-	  Dc.M(src,r);
-	}
-	FGrid->Barrier();
-	double t1=usecond();
-	uint64_t ncall = 500;
-
-	FGrid->Broadcast(0,&ncall,sizeof(ncall));
-
-	//	std::cout << GridLogMessage << " Estimate " << ncall << " calls per second"<<std::endl;
-
-	time_statistics timestat;
-	std::vector<double> t_time(ncall);
-	for(uint64_t i=0;i<ncall;i++){
-	  t0=usecond();
-	  Dc.M(src,r);
-	  t1=usecond();
-	  t_time[i] = t1-t0;
-	}
-	FGrid->Barrier();
-	
-	double volume=1;  for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
-	double flops=(1344+ 24+6*6*8*2)*volume;
-	double mf_hi, mf_lo, mf_err;
-	
-	timestat.statistics(t_time);
-	mf_hi = flops/timestat.min;
-	mf_lo = flops/timestat.max;
-	mf_err= flops/timestat.min * timestat.err/timestat.mean;
-
-	mflops = flops/timestat.mean;
-	mflops_all.push_back(mflops);
-	if ( mflops_best == 0   ) mflops_best = mflops;
-	if ( mflops_worst== 0   ) mflops_worst= mflops;
-	if ( mflops>mflops_best ) mflops_best = mflops;
-	if ( mflops<mflops_worst) mflops_worst= mflops;
-	
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Dclov mflop/s =   "<< mflops << " ("<<mf_err<<") " << mf_lo<<"-"<<mf_hi <<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Dclov mflop/s per rank   "<< mflops/NP<<std::endl;
-	std::cout<<GridLogMessage << std::fixed << std::setprecision(1)<<"Dclov mflop/s per node   "<< mflops/NN<<std::endl;
-      
-      }
-
-      std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Best  mflop/s        =   "<< mflops_best << " ; " << mflops_best/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage << L<<"^4  Deo Worst mflop/s        =   "<< mflops_worst<< " ; " << mflops_worst/NN<<" per node " <<std::endl;
-      std::cout<<GridLogMessage <<fmt << std::endl;
-      std::cout<<GridLogMessage ;
-
-      for(int i=0;i<mflops_all.size();i++){
-	std::cout<<mflops_all[i]/NN<<" ; " ;
-      }
-      std::cout<<std::endl;
-    }
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    return mflops_best;
-  }
-};
-
-
-
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  if (GlobalSharedMemory::WorldRank==0) { 
-    FP = fopen("Benchmark_usqcd.csv","w");
-  } else {
-    FP = fopen("/dev/null","w");
-  }
-
-  CartesianCommunicator::SetCommunicatorPolicy(CartesianCommunicator::CommunicatorPolicySequential);
-  LebesgueOrder::Block = std::vector<int>({2,2,2,2});
-
-  Benchmark::Decomposition();
-
-  int do_su4=0;
-  int do_memory=1;
-  int do_comms =1;
-  int do_blas  =1;
-
-  int sel=4;
-  std::vector<int> L_list({8,12,16,24,32});
-  int selm1=sel-1;
-
-  std::vector<double> clover;
-  std::vector<double> dwf4;
-  std::vector<double> staggered;
-
-  int Ls=1;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << " Clover dslash 4D vectorised (temporarily Wilson)" <<std::endl;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  for(int l=0;l<L_list.size();l++){
-    clover.push_back(Benchmark::DWF(1,L_list[l]));
-  }
-
-  Ls=12;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << " Domain wall dslash 4D vectorised" <<std::endl;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  for(int l=0;l<L_list.size();l++){
-    double result = Benchmark::DWF(Ls,L_list[l]) ;
-    dwf4.push_back(result);
-  }
-
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << " Improved Staggered dslash 4D vectorised" <<std::endl;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  for(int l=0;l<L_list.size();l++){
-    double result = Benchmark::Staggered(L_list[l]) ;
-    staggered.push_back(result);
-  }
-
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << " Summary table Ls="<<Ls <<std::endl;
-  std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-  std::cout<<GridLogMessage << "L \t\t Clover \t\t DWF4 \t\t Staggered" <<std::endl;
-  for(int l=0;l<L_list.size();l++){
-    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_memory ) {
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Memory benchmark " <<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    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;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    Benchmark::SU4();
-  }
-  
-  if ( do_comms ) {
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    std::cout<<GridLogMessage << " Communications benchmark " <<std::endl;
-    std::cout<<GridLogMessage << "=================================================================================="<<std::endl;
-    Benchmark::Comms();
-  }
-
-    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;
-
-  Grid_finalize();
-  fclose(FP);
-}

bootstrap.sh

@@ -1,12 +1,12 @@
 #!/usr/bin/env bash
 set -e
 
-EIGEN_URL='https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.bz2'
-EIGEN_SHA256SUM='b4c198460eba6f28d34894e3a5710998818515104d6e74e5cc331ce31e46e626'
+EIGEN_URL='https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2'
+EIGEN_SHA256SUM='685adf14bd8e9c015b78097c1dc22f2f01343756f196acdc76a678e1ae352e11'
 
 
 echo "-- deploying Eigen source..."
-ARC=$(basename ${EIGEN_URL})
+ARC=`basename ${EIGEN_URL}`
 wget ${EIGEN_URL} --no-check-certificate
 if command -v sha256sum; then
    echo "$EIGEN_SHA256SUM  $(basename "$EIGEN_URL")" \
@@ -14,8 +14,13 @@ if command -v sha256sum; then
 else
    echo "WARNING: could not verify checksum, please install sha256sum" >&2
 fi
-./scripts/update_eigen.sh "${ARC}"
-rm "${ARC}"
+./scripts/update_eigen.sh ${ARC}
+rm ${ARC}
+# patch for non-portable includes in Eigen 3.3.5
+# apparently already fixed in Eigen HEAD so it should not be 
+# a problem in the future (A.P.)
+patch Eigen/unsupported/Eigen/CXX11/Tensor scripts/eigen-3.3.5.Tensor.patch
+
 echo '-- generating Make.inc files...'
 ./scripts/filelist
 echo '-- generating configure script...'

configure.ac

@@ -226,14 +226,23 @@ case ${ac_SFW_FP16} in
 esac
 
 ############### Default to accelerator cshift, but revert to host if UCX is buggy or other reasons
-AC_ARG_ENABLE([accelerator-aware-mpi],
-    [AS_HELP_STRING([--enable-accelerator-aware-mpi=yes|no],[run mpi transfers from device])],
-    [ac_ACCELERATOR_AWARE_MPI=${enable_accelerator_aware_mpi}], [ac_ACCELERATOR_AWARE_MPI=yes])
+AC_ARG_ENABLE([accelerator-cshift],
+    [AS_HELP_STRING([--enable-accelerator-cshift=yes|no],[run cshift on the device])],
+    [ac_ACC_CSHIFT=${enable_accelerator_cshift}], [ac_ACC_CSHIFT=yes])
 
-case ${ac_ACCELERATOR_AWARE_MPI} in
+AC_ARG_ENABLE([ucx-buggy],
+    [AS_HELP_STRING([--enable-ucx-buggy=yes|no],[enable workaround for UCX device buffer bugs])],
+    [ac_UCXBUGGY=${enable_ucx_buggy}], [ac_UCXBUGGY=no])
+
+case ${ac_UCXBUGGY} in
     yes)
-      AC_DEFINE([ACCELERATOR_CSHIFT],[1],[ Cshift runs on host])
-      AC_DEFINE([ACCELERATOR_AWARE_MPI],[1],[ Stencil can use device pointers]);;
+    ac_ACC_CSHIFT=no;;
+    *);;
+esac
+
+case ${ac_ACC_CSHIFT} in
+    yes)
+      AC_DEFINE([ACCELERATOR_CSHIFT],[1],[ UCX device buffer bugs are not present]);;
     *);;
 esac
 

scripts/eigen-3.3.5.Tensor.patch

@@ -0,0 +1,19 @@
+--- ./Eigen/unsupported/Eigen/CXX11/Tensor	2018-07-23 10:33:42.000000000 +0100
++++ Tensor	2018-08-28 16:15:56.000000000 +0100
+@@ -25,7 +25,7 @@
+ #include <utility>
+ #endif
+ 
+-#include <Eigen/src/Core/util/DisableStupidWarnings.h>
++#include "../../../Eigen/src/Core/util/DisableStupidWarnings.h"
+ 
+ #include "../SpecialFunctions"
+ #include "src/util/CXX11Meta.h"
+@@ -147,6 +147,6 @@
+ 
+ #include "src/Tensor/TensorIO.h"
+ 
+-#include <Eigen/src/Core/util/ReenableStupidWarnings.h>
++#include "../../../Eigen/src/Core/util/ReenableStupidWarnings.h"
+ 
+ //#endif // EIGEN_CXX11_TENSOR_MODULE

systems/Aurora/benchmarks/bench1024.pbs

@@ -25,16 +25,12 @@ export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
 
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-#export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
+export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
 export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
 export MPICH_OFI_NIC_POLICY=GPU
-export FI_CXI_CQ_FILL_PERCENT=10
-export FI_CXI_DEFAULT_CQ_SIZE=262144
-#export FI_CXI_DEFAULT_CQ_SIZE=131072
-#export FI_CXI_CQ_FILL_PERCENT=20
 
 # 12 ppn, 32 nodes, 384 ranks
 #
@@ -49,12 +45,12 @@ CMD="mpiexec -np 12288 -ppn 12  -envall \
 	     ./gpu_tile_compact.sh \
 	     ./Benchmark_dwf_fp32 --mpi 8.8.8.24 --grid 128.128.128.384 \
 		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32 --comms-overlap"
-$CMD | tee 1024node.dwf.small.cq
+$CMD | tee 1024node.dwf.small
 
 CMD="mpiexec -np 12288 -ppn 12  -envall \
 	     ./gpu_tile_compact.sh \
 	     ./Benchmark_dwf_fp32 --mpi 16.8.8.12 --grid 256.256.256.384 \
 		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32 --comms-overlap"
-$CMD | tee 1024node.dwf.cq
+$CMD | tee 1024node.dwf
 
 

systems/Aurora/benchmarks/bench12.pbs

@@ -17,7 +17,6 @@ source ../sourceme.sh
 export OMP_NUM_THREADS=3
 export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
 
-
 #unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE
 #unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE
 #unset MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST
@@ -36,25 +35,11 @@ CMD="mpiexec -np 24 -ppn 12  -envall \
 	     ./Benchmark_comms_host_device --mpi 2.3.2.2 --grid 32.24.32.192 \
 		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32"
 
-#$CMD 
+$CMD 
 
 CMD="mpiexec -np 24 -ppn 12  -envall \
 	     ./gpu_tile_compact.sh \
 	     ./Benchmark_dwf_fp32 --mpi 2.3.2.2 --grid 64.96.64.64 --comms-overlap \
 		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32"
 
-#$CMD 
-
-CMD="mpiexec -np 1 -ppn 1  -envall \
-	     ./gpu_tile_compact.sh \
-	     ./Benchmark_dwf --mpi 1.1.1.1 --grid 16.32.32.32 --comms-sequential \
-		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32"
-
-$CMD 
-
-CMD="mpiexec -np 1 -ppn 1  -envall \
-	     ./gpu_tile_compact.sh \
-	     ./Benchmark_dwf_fp32 --mpi 1.1.1.1 --grid 16.32.32.32 --comms-sequential \
-		--shm-mpi 1 --shm 2048 --device-mem 32000 --accelerator-threads 32"
-
 $CMD 

systems/Aurora/config-command

@@ -1,16 +1,16 @@
-
+TOOLS=$HOME/tools
 ../../configure \
 	--enable-simd=GPU \
 	--enable-gen-simd-width=64 \
 	--enable-comms=mpi-auto \
+	--enable-accelerator-cshift \
 	--disable-gparity \
 	--disable-fermion-reps \
 	--enable-shm=nvlink \
 	--enable-accelerator=sycl \
-	--enable-accelerator-aware-mpi=no\
 	--enable-unified=no \
 	MPICXX=mpicxx \
 	CXX=icpx \
-	LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel -lsycl" \
-	CXXFLAGS="-fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -qmkl=parallel"
+	LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L$TOOLS/lib64/" \
+	CXXFLAGS="-fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -I$TOOLS/include"
 

systems/Aurora/sourceme-sunspot-deterministic.sh

@@ -1,2 +0,0 @@
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu

systems/Aurora/sourceme.sh

@@ -3,24 +3,10 @@
 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
-# -ftarget-register-alloc-mode=pvc:auto
-#
-
 export HTTP_PROXY=http://proxy.alcf.anl.gov:3128
 export HTTPS_PROXY=http://proxy.alcf.anl.gov:3128
 export http_proxy=http://proxy.alcf.anl.gov:3128
 export https_proxy=http://proxy.alcf.anl.gov:3128
 #export MPIR_CVAR_CH4_OFI_ENABLE_HMEM=1
 git config --global http.proxy http://proxy.alcf.anl.gov:3128
-
-export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
+ 

systems/Aurora/tests/repro128.pbs

@@ -1,41 +0,0 @@
-#!/bin/bash
-
-## qsub -q EarlyAppAccess -A Aurora_Deployment -I -l select=1 -l walltime=60:00
-
-#PBS -q EarlyAppAccess
-#PBS -l select=128
-#PBS -l walltime=02:00:00
-#PBS -A LatticeQCD_aesp_CNDA
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-cd $PBS_O_WORKDIR
-
-source ../sourceme.sh
-
-cat $PBS_NODEFILE
-
-export OMP_NUM_THREADS=3
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE
-#unset MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST
-
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-export MPICH_OFI_NIC_POLICY=GPU
-
-# 12 ppn, 16 nodes, 192 ranks
-# 12 ppn, 128 nodes, 1536 ranks
-CMD="mpiexec -np 1536 -ppn 12  -envall \
-	     ./gpu_tile_compact.sh \
-	     ./Test_dwf_mixedcg_prec --mpi 4.4.4.24 --grid 128.128.128.384 \
-		--shm-mpi 1 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 7000 --comms-overlap "
-$CMD 

systems/Aurora/tests/repro16.pbs

@@ -1,61 +0,0 @@
-#!/bin/bash
-
-## qsub -q EarlyAppAccess -A Aurora_Deployment -I -l select=1 -l walltime=60:00
-
-#PBS -l select=16:system=sunspot,place=scatter
-#PBS -A LatticeQCD_aesp_CNDA
-#PBS -l walltime=01:00:00
-#PBS -N dwf
-#PBS -k doe
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-cd $PBS_O_WORKDIR
-
-#source ../sourceme.sh
-
-cat $PBS_NODEFILE
-
-#export MPICH_COLL_SYNC=1
-#export MPICH_ENV_DISPLAY=1
-export MPICH_
-export OMP_NUM_THREADS=3
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu
-#export LD_LIBRARY_PATH=/soft/restricted/CNDA/updates/2023.05.15.001/oneapi/compiler/eng-20230512/compiler/linux/lib/:$LD_LIBRARY_PATH
-
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE
-#unset MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST
-export MPIR_CVAR_ALLREDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_REDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_ALLREDUCE_INTRA_ALGORITHM=recursive_doubling
-unset MPIR_CVAR_CH4_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_CH4_POSIX_COLL_SELECTION_TUNING_JSON_FILE
-
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-export MPICH_OFI_NIC_POLICY=GPU
-
-DIR=repro.$PBS_JOBID
-mkdir $DIR
-cd $DIR
-
-CMD="mpiexec -np 192 -ppn 12  -envall \
-	     ../gpu_tile_compact.sh \
-	     ../Test_dwf_mixedcg_prec --mpi 2.4.4.6 --grid 64.128.128.192 \
-		--shm-mpi 1 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 3000 --debug-stdout --log Message,Iterative"
-#--comms-overlap
-$CMD 
-
-grep Oops Grid.stderr.* > failures.$PBS_JOBID
-rm core.*
-

systems/Aurora/tests/repro1gpu.pbs

@@ -1,82 +0,0 @@
-#!/bin/bash
-
-#PBS -l select=16:system=sunspot,place=scatter
-#PBS -A LatticeQCD_aesp_CNDA
-#PBS -l walltime=02:00:00
-#PBS -N repro1gpu
-#PBS -k doe
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu
-
-# 56 cores / 6 threads ~9
-export OMP_NUM_THREADS=6
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-export MPICH_OFI_NIC_POLICY=GPU
-
-export MPIR_CVAR_ALLREDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_REDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_ALLREDUCE_INTRA_ALGORITHM=recursive_doubling
-unset MPIR_CVAR_CH4_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_CH4_POSIX_COLL_SELECTION_TUNING_JSON_FILE
-export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
-
-cd $PBS_O_WORKDIR
-
-NN=`cat $PBS_NODEFILE | wc -l`
-echo $PBS_NODEFILE
-cat $PBS_NODEFILE
-
-echo $NN nodes in node file
-for n in `eval echo {1..$NN}`
-do
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-echo Node $n is $THIS_NODE
-
-
-for g in {0..11}
-do
-export NUMA_MAP=(0 0 0 1 1 1 0 0 0 1 1 1 )
-export TILE_MAP=(0 0 0 0 0 0 1 1 1 1 1 1 )
-export  GPU_MAP=(0 1 2 3 4 5 0 1 2 3 4 5 )
-
-export numa=${NUMA_MAP[$g]}
-export gpu_id=${GPU_MAP[$g]}
-export tile_id=${TILE_MAP[$g]}
-export gpu=$gpu_id.$tile_id
-
-cd $PBS_O_WORKDIR
-
-DIR=repro.1gpu.$PBS_JOBID/node-$n-$THIS_NODE-GPU-$gpu
-mkdir -p $DIR
-cd $DIR
-
-echo $THIS_NODE > nodefile
-echo $gpu > gpu
-
-export ZE_AFFINITY_MASK=$gpu
-export ONEAPI_DEVICE_FILTER=gpu,level_zero
-
-CMD="mpiexec -np 1 -ppn 1  -envall --hostfile nodefile \
-	     numactl -N $numa -m $numa ../../Test_dwf_mixedcg_prec --mpi 1.1.1.1 --grid 16.16.32.32 \
-		--shm-mpi 0 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message"
-echo $CMD
-$CMD &
-
-done
-done
-
-wait
-

systems/Aurora/tests/reproN.pbs

@@ -1,98 +0,0 @@
-#!/bin/bash
-
-#PBS -l select=32:system=sunspot,place=scatter
-#PBS -A LatticeQCD_aesp_CNDA
-#PBS -l walltime=02:00:00
-#PBS -N reproN
-#PBS -k doe
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu
-
-# 56 cores / 6 threads ~9
-export OMP_NUM_THREADS=6
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-#export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-
-export SYCL_PI_LEVEL_ZERO_USE_IMMEDIATE_COMMANDLISTS=1
-export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE=1
-export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE_FOR_D2D_COPY=1
-export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"
-
-export GRID_PRINT_ENTIRE_LOG=0
-export GRID_CHECKSUM_RECV_BUF=0
-export GRID_CHECKSUM_SEND_BUF=0
-
-export MPICH_OFI_NIC_POLICY=GPU
-
-export MPIR_CVAR_ALLREDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_REDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_ALLREDUCE_INTRA_ALGORITHM=recursive_doubling
-unset MPIR_CVAR_CH4_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_CH4_POSIX_COLL_SELECTION_TUNING_JSON_FILE
-
-cd $PBS_O_WORKDIR
-
-NN=`cat $PBS_NODEFILE | wc -l`
-echo $PBS_NODEFILE
-cat $PBS_NODEFILE
-
-echo $NN nodes in node file
-for n in `eval echo {1..$NN}`
-do
-
-cd $PBS_O_WORKDIR
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-echo Node $n is $THIS_NODE
-
-DIR=reproN.$PBS_JOBID/node-$n-$THIS_NODE
-
-mkdir -p $DIR
-cd $DIR
-
-echo $THIS_NODE > nodefile
-
-#CMD="mpiexec -np 12 -ppn 12  -envall --hostfile nodefile \
-#	     ../../gpu_tile_compact.sh \
-#	     ../../Test_dwf_mixedcg_prec --mpi 1.2.2.3 --grid 32.64.64.96 \
-#		--shm-mpi 0 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message --comms-overlap"
-
-CMD="mpiexec -np 12 -ppn 12  -envall --hostfile nodefile \
-	     ../../gpu_tile_compact.sh \
-	     ../../Test_dwf_mixedcg_prec --mpi 1.2.2.3 --grid 32.64.64.96 \
-		--shm-mpi 1 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message --comms-overlap"
-
-echo $CMD > command-line
-env > environment
-$CMD &
-
-done
-
-# Suspicious wait is allowing jobs to collide and knock out
-#wait
-
-sleep 6500
-
-for n in ` eval echo {1..$NN} `
-do
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-DIR=reproN.$PBS_JOBID/node-$n-$THIS_NODE
-
-cd $DIR
-
-grep Oops Grid.stderr.* > failures.$PBS_JOBID
-rm core.*
-
-done

systems/Aurora/tests/solver/stag16.pbs

@@ -1,40 +0,0 @@
-#!/bin/bash
-
-## qsub -q EarlyAppAccess -A Aurora_Deployment -I -l select=1 -l walltime=60:00
-
-#PBS -q EarlyAppAccess
-#PBS -l select=16
-#PBS -l walltime=01:00:00
-#PBS -A LatticeQCD_aesp_CNDA
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-cd $PBS_O_WORKDIR
-
-source ../../sourceme.sh
-
-cat $PBS_NODEFILE
-
-export OMP_NUM_THREADS=3
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE
-#unset MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE
-#unset MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST
-
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-export MPICH_OFI_NIC_POLICY=GPU
-
-# 12 ppn, 16 nodes, 192 ranks
-CMD="mpiexec -np 192 -ppn 12  -envall \
-	     ./gpu_tile_compact.sh \
-	     ./Test_staggered_cg_prec --mpi 2.4.4.6 --grid 128.128.128.192 \
-	     --shm-mpi 1 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 3000 --comms-overlap"
-$CMD 

systems/Booster/benchmarks/Benchmark_usqcd.csv

@@ -1,70 +0,0 @@
-Memory Bandwidth
-
-Bytes, GB/s per node
-3145728, 225.900365
-50331648, 2858.859504
-254803968, 4145.556367
-805306368, 4905.772480
-1966080000, 4978.312557
-
-
-GEMM
-
- M, N, K, BATCH, GF/s per rank
-16, 8, 16, 256, 1.713639
-16, 16, 16, 256, 288.268316
-16, 32, 16, 256, 597.053950
-32, 8, 32, 256, 557.382591
-32, 16, 32, 256, 1100.145311
-32, 32, 32, 256, 1885.080449
-64, 8, 64, 256, 1725.163599
-64, 16, 64, 256, 3389.336566
-64, 32, 64, 256, 4168.252422
-16, 8, 256, 256, 1326.262134
-16, 16, 256, 256, 2318.095475
-16, 32, 256, 256, 3555.436503
-32, 8, 256, 256, 1920.139170
-32, 16, 256, 256, 3486.174753
-32, 32, 256, 256, 5320.821724
-64, 8, 256, 256, 2539.597502
-64, 16, 256, 256, 5003.456775
-64, 32, 256, 256, 7837.531562
-8, 256, 16, 256, 1427.848170
-16, 256, 16, 256, 2222.147815
-32, 256, 16, 256, 2877.121715
-8, 256, 32, 256, 1922.890086
-16, 256, 32, 256, 3199.469082
-32, 256, 32, 256, 4845.405343
-8, 256, 64, 256, 2639.483343
-16, 256, 64, 256, 5012.800299
-32, 256, 64, 256, 7216.006882
-
-
-
-Communications
-
-Packet bytes, direction, GB/s per node
-4718592, 2, 206.570734
-4718592, 3, 207.501847
-4718592, 6, 189.730277
-4718592, 7, 204.301218
-15925248, 2, 307.882997
-15925248, 3, 287.901076
-15925248, 6, 295.603109
-15925248, 7, 300.682033
-37748736, 2, 331.740364
-37748736, 3, 338.610627
-37748736, 6, 332.580657
-37748736, 7, 336.336579
-
-
-Per node summary table
-
-L , Wilson, DWF4, Staggered, GF/s per node
-
-8 , 16, 1165, 10
-12 , 473, 4901, 163
-16 , 1436, 8464, 442
-24 , 4133, 10139, 1530
-32 , 5726, 11487, 2518
-

systems/Booster/config-command

@@ -5,12 +5,10 @@ LIME=/p/home/jusers/boyle2/juwels/gm2dwf/boyle/
     --enable-gen-simd-width=64 \
     --enable-shm=nvlink \
     --enable-accelerator=cuda \
-    --disable-gparity \
-    --disable-fermion-reps \
     --with-lime=$LIME \
-    --enable-accelerator-cshift \
+    --disable-accelerator-cshift \
     --disable-unified \
     CXX=nvcc \
     LDFLAGS="-cudart shared " \
-    CXXFLAGS="-ccbin mpicxx -gencode arch=compute_80,code=sm_80 -std=c++17 -cudart shared -lcublas"
+    CXXFLAGS="-ccbin mpicxx -gencode arch=compute_80,code=sm_80 -std=c++14 -cudart shared"
 

systems/Booster/sourceme.sh

@@ -1,5 +1,5 @@
-module load GCC
-module load GMP
-module load MPFR
-module load OpenMPI
-module load CUDA
+module load GCC/9.3.0       
+module load  GMP/6.2.0   
+module load MPFR/4.1.0     
+module load OpenMPI/4.1.0rc1  
+module load CUDA/11.3

systems/Frontier/config-command

@@ -16,7 +16,7 @@ CLIME=`spack find --paths c-lime@2-3-9 | grep c-lime| cut -c 15-`
 --disable-fermion-reps \
 CXX=hipcc MPICXX=mpicxx \
 CXXFLAGS="-fPIC -I{$ROCM_PATH}/include/ -I${MPICH_DIR}/include -L/lib64 -fgpu-sanitize" \
- LDFLAGS="-L/lib64 -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64  -lhipblas -lrocblas"
+ LDFLAGS="-L/lib64 -L${MPICH_DIR}/lib -lmpi -L${CRAY_MPICH_ROOTDIR}/gtl/lib -lmpi_gtl_hsa -lamdhip64 "
 
 
 

systems/PVC-OEM/setup.sh

@@ -1,5 +1,3 @@
 export https_proxy=http://proxy-chain.intel.com:911
 module load intel-release
 module load intel/mpich
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-export SYCL_PROGRAM_COMPILE_OPTIONS="-ze-opt-large-register-file"

systems/Sunspot/config-command

@@ -1,4 +1,4 @@
-
+TOOLS=$HOME/tools
 ../../configure \
 	--enable-simd=GPU \
 	--enable-gen-simd-width=64 \
@@ -11,6 +11,6 @@
 	--enable-unified=no \
 	MPICXX=mpicxx \
 	CXX=icpx \
-	LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L${MKLROOT}/lib -qmkl=parallel -lsycl" \
-	CXXFLAGS="-fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -qmkl=parallel"
+	LDFLAGS="-fiopenmp -fsycl -fsycl-device-code-split=per_kernel -fsycl-device-lib=all -lze_loader -L$TOOLS/lib64/" \
+	CXXFLAGS="-fiopenmp -fsycl-unnamed-lambda -fsycl -I$INSTALL/include -Wno-tautological-compare -I$HOME/ -I$TOOLS/include"
 

systems/Sunspot/sourceme.sh

@@ -1,2 +0,0 @@
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu

systems/Sunspot/tests/repro1gpu.pbs

@@ -1,81 +0,0 @@
-#!/bin/bash
-
-#PBS -l select=16:system=sunspot,place=scatter
-#PBS -A LatticeQCD_aesp_CNDA
-#PBS -l walltime=02:00:00
-#PBS -N repro1gpu
-#PBS -k doe
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu
-
-# 56 cores / 6 threads ~9
-export OMP_NUM_THREADS=6
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-export MPICH_OFI_NIC_POLICY=GPU
-
-export MPIR_CVAR_ALLREDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_REDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_ALLREDUCE_INTRA_ALGORITHM=recursive_doubling
-unset MPIR_CVAR_CH4_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_CH4_POSIX_COLL_SELECTION_TUNING_JSON_FILE
-
-cd $PBS_O_WORKDIR
-
-NN=`cat $PBS_NODEFILE | wc -l`
-echo $PBS_NODEFILE
-cat $PBS_NODEFILE
-
-echo $NN nodes in node file
-for n in `eval echo {1..$NN}`
-do
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-echo Node $n is $THIS_NODE
-
-
-for g in {0..11}
-do
-export NUMA_MAP=(0 0 0 1 1 1 0 0 0 1 1 1 )
-export TILE_MAP=(0 0 0 0 0 0 1 1 1 1 1 1 )
-export  GPU_MAP=(0 1 2 3 4 5 0 1 2 3 4 5 )
-
-export numa=${NUMA_MAP[$g]}
-export gpu_id=${GPU_MAP[$g]}
-export tile_id=${TILE_MAP[$g]}
-export gpu=$gpu_id.$tile_id
-
-cd $PBS_O_WORKDIR
-
-DIR=repro.1gpu.$PBS_JOBID/node-$n-$THIS_NODE-GPU-$gpu
-mkdir -p $DIR
-cd $DIR
-
-echo $THIS_NODE > nodefile
-echo $gpu > gpu
-
-export ZE_AFFINITY_MASK=$gpu
-export ONEAPI_DEVICE_FILTER=gpu,level_zero
-
-CMD="mpiexec -np 1 -ppn 1  -envall --hostfile nodefile \
-	     numactl -N $numa -m $numa ../../Test_dwf_mixedcg_prec --mpi 1.1.1.1 --grid 16.16.32.32 \
-		--shm-mpi 0 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message"
-echo $CMD
-$CMD &
-
-done
-done
-
-wait
-

systems/Sunspot/tests/reproN.pbs

@@ -1,97 +0,0 @@
-#!/bin/bash
-
-#PBS -l select=32:system=sunspot,place=scatter
-#PBS -A LatticeQCD_aesp_CNDA
-#PBS -l walltime=02:00:00
-#PBS -N reproN
-#PBS -k doe
-
-#export OMP_PROC_BIND=spread
-#unset OMP_PLACES
-
-module load oneapi/eng-compiler/2023.05.15.003
-module load mpich/51.2/icc-all-deterministic-pmix-gpu
-
-# 56 cores / 6 threads ~9
-export OMP_NUM_THREADS=6
-export MPIR_CVAR_CH4_OFI_ENABLE_GPU_PIPELINE=1
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_D2H_ENGINE_TYPE=0
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_H2D_ENGINE_TYPE=0
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_BUFFER_SZ=1048576
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_THRESHOLD=131072
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_NUM_BUFFERS_PER_CHUNK=16
-#export MPIR_CVAR_CH4_OFI_GPU_PIPELINE_MAX_NUM_BUFFERS=16
-#export MPIR_CVAR_GPU_USE_IMMEDIATE_COMMAND_LIST=1
-
-export SYCL_PI_LEVEL_ZERO_USE_IMMEDIATE_COMMANDLISTS=1
-export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE=1
-export SYCL_PI_LEVEL_ZERO_USE_COPY_ENGINE_FOR_D2D_COPY=1
-
-export GRID_PRINT_ENTIRE_LOG=0
-export GRID_CHECKSUM_RECV_BUF=1
-export GRID_CHECKSUM_SEND_BUF=0
-
-export MPICH_OFI_NIC_POLICY=GPU
-
-export MPIR_CVAR_ALLREDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_REDUCE_DEVICE_COLLECTIVE=0
-export MPIR_CVAR_ALLREDUCE_INTRA_ALGORITHM=recursive_doubling
-unset MPIR_CVAR_CH4_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_COLL_SELECTION_TUNING_JSON_FILE
-unset MPIR_CVAR_CH4_POSIX_COLL_SELECTION_TUNING_JSON_FILE
-
-cd $PBS_O_WORKDIR
-
-NN=`cat $PBS_NODEFILE | wc -l`
-echo $PBS_NODEFILE
-cat $PBS_NODEFILE
-
-echo $NN nodes in node file
-for n in `eval echo {1..$NN}`
-do
-
-cd $PBS_O_WORKDIR
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-echo Node $n is $THIS_NODE
-
-DIR=reproN.$PBS_JOBID/node-$n-$THIS_NODE
-
-mkdir -p $DIR
-cd $DIR
-
-echo $THIS_NODE > nodefile
-
-#CMD="mpiexec -np 12 -ppn 12  -envall --hostfile nodefile \
-#	     ../../gpu_tile_compact.sh \
-#	     ../../Test_dwf_mixedcg_prec --mpi 1.2.2.3 --grid 32.64.64.96 \
-#		--shm-mpi 0 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message --comms-overlap"
-
-CMD="mpiexec -np 12 -ppn 12  -envall --hostfile nodefile \
-	     ../../gpu_tile_compact.sh \
-	     ../../Test_dwf_mixedcg_prec --mpi 1.2.2.3 --grid 32.64.64.96 \
-		--shm-mpi 1 --shm 4096 --device-mem 32000 --accelerator-threads 32 --seconds 6000 --debug-stdout --log Message --comms-overlap"
-
-echo $CMD > command-line
-env > environment
-$CMD &
-
-done
-
-# Suspicious wait is allowing jobs to collide and knock out
-#wait
-
-sleep 6500
-
-for n in ` eval echo {1..$NN} `
-do
-
-THIS_NODE=`head -n$n $PBS_NODEFILE | tail -n1 `
-DIR=reproN.$PBS_JOBID/node-$n-$THIS_NODE
-
-cd $DIR
-
-grep Oops Grid.stderr.* > failures.$PBS_JOBID
-rm core.*
-
-done

systems/mac-arm/config-command-mpi

@@ -1,3 +1,4 @@
-CXXFLAGS=-I/opt/local/include LDFLAGS=-L/opt/local/lib/ CXX=c++-13 MPICXX=mpicxx ../../configure --enable-simd=GEN --enable-comms=mpi-auto --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-fermion-reps --disable-gparity --disable-debug 
+BREW=/opt/local/
+MPICXX=mpicxx ../../configure --enable-simd=GEN --enable-comms=mpi-auto --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-fermion-reps --disable-gparity --disable-debug
 
 

tests/Test_dwf_mixedcg_prec.cc

@@ -30,60 +30,27 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 using namespace std;
 using namespace Grid;
 
-#ifndef HOST_NAME_MAX
-#define HOST_NAME_MAX _POSIX_HOST_NAME_MAX
-#endif
-
-
-NAMESPACE_BEGIN(Grid);
-template<class Matrix,class Field>
-  class SchurDiagMooeeOperatorParanoid :  public SchurOperatorBase<Field> {
- public:
-    Matrix &_Mat;
-    SchurDiagMooeeOperatorParanoid (Matrix &Mat): _Mat(Mat){};
-    virtual  void Mpc      (const Field &in, Field &out) {
-      Field tmp(in.Grid());
-      tmp.Checkerboard() = !in.Checkerboard();
-      //      std::cout <<" Mpc starting"<<std::endl;
-
-      RealD nn = norm2(in); // std::cout <<" Mpc Prior to dslash norm is "<<nn<<std::endl;
-      _Mat.Meooe(in,tmp);
-      nn = norm2(tmp); //std::cout <<" Mpc Prior to Mooeinv "<<nn<<std::endl;
-      _Mat.MooeeInv(tmp,out);
-      nn = norm2(out); //std::cout <<" Mpc Prior to dslash norm is "<<nn<<std::endl;
-      _Mat.Meooe(out,tmp);
-      nn = norm2(tmp); //std::cout <<" Mpc Prior to Mooee "<<nn<<std::endl;
-      _Mat.Mooee(in,out);
-      nn = norm2(out); //std::cout <<" Mpc Prior to axpy "<<nn<<std::endl;
-      axpy(out,-1.0,tmp,out);
-    }
-    virtual void MpcDag   (const Field &in, Field &out){
-      Field tmp(in.Grid());
-      //      std::cout <<" MpcDag starting"<<std::endl;
-      RealD nn = norm2(in);// std::cout <<" MpcDag Prior to dslash norm is "<<nn<<std::endl;
-      _Mat.MeooeDag(in,tmp);
-      _Mat.MooeeInvDag(tmp,out);
-      nn = norm2(out);// std::cout <<" MpcDag Prior to dslash norm is "<<nn<<std::endl;
-      _Mat.MeooeDag(out,tmp);
-      nn = norm2(tmp);// std::cout <<" MpcDag Prior to Mooee "<<nn<<std::endl;
-      _Mat.MooeeDag(in,out);
-      nn = norm2(out);// std::cout <<" MpcDag Prior to axpy "<<nn<<std::endl;
-      axpy(out,-1.0,tmp,out);
-    }
+template<class d>
+struct scal {
+  d internal;
 };
 
-NAMESPACE_END(Grid);
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
 
 int main (int argc, char ** argv)
 {
-  char hostname[HOST_NAME_MAX+1];
-  gethostname(hostname, HOST_NAME_MAX+1);
-  std::string host(hostname);
-  
   Grid_init(&argc,&argv);
 
   const int Ls=12;
 
+  std::cout << GridLogMessage << "::::: NB: to enable a quick bit reproducibility check use the --checksums flag. " << std::endl;
+
+  { 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
@@ -122,17 +89,10 @@ int main (int argc, char ** argv)
   result_o_2.Checkerboard() = Odd;
   result_o_2 = Zero();
 
-  SchurDiagMooeeOperatorParanoid<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
-  SchurDiagMooeeOperatorParanoid<DomainWallFermionF,LatticeFermionF> HermOpEO_f(Ddwf_f);
-
-  int nsecs=600;
-  if( GridCmdOptionExists(argv,argv+argc,"--seconds") ){
-    std::string arg = GridCmdOptionPayload(argv,argv+argc,"--seconds");
-    GridCmdOptionInt(arg,nsecs);
-  }
-  
-  std::cout << GridLogMessage << "::::::::::::: Starting mixed CG for "<<nsecs <<" seconds" << std::endl;
+  SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermionD> HermOpEO(Ddwf);
+  SchurDiagMooeeOperator<DomainWallFermionF,LatticeFermionF> HermOpEO_f(Ddwf_f);
 
+  std::cout << GridLogMessage << "::::::::::::: Starting mixed CG" << std::endl;
   MixedPrecisionConjugateGradient<LatticeFermionD,LatticeFermionF> mCG(1.0e-8, 10000, 50, FrbGrid_f, HermOpEO_f, HermOpEO);
   double t1,t2,flops;
   double MdagMsiteflops = 1452; // Mobius (real coeffs)
@@ -141,26 +101,7 @@ int main (int argc, char ** argv)
   std:: cout << " MdagM site flops = "<< 4*MdagMsiteflops<<std::endl;
   std:: cout << " CG    site flops = "<< CGsiteflops <<std::endl;
   int iters;
-
-  time_t start = time(NULL);
-
-  FlightRecorder::ContinueOnFail = 0;
-  FlightRecorder::PrintEntireLog = 0;
-  FlightRecorder::ChecksumComms  = 1;
-  FlightRecorder::ChecksumCommsSend=0;
-
-  if(char *s=getenv("GRID_PRINT_ENTIRE_LOG"))  FlightRecorder::PrintEntireLog     = atoi(s);
-  if(char *s=getenv("GRID_CHECKSUM_RECV_BUF")) FlightRecorder::ChecksumComms      = atoi(s);
-  if(char *s=getenv("GRID_CHECKSUM_SEND_BUF")) FlightRecorder::ChecksumCommsSend  = atoi(s);
-
-  int iter=0;
-  do {
-    if ( iter == 0 ) {
-      FlightRecorder::SetLoggingMode(FlightRecorder::LoggingModeRecord);
-    } else {
-      FlightRecorder::SetLoggingMode(FlightRecorder::LoggingModeVerify);
-    }
-    std::cerr << "******************* SINGLE PRECISION SOLVE "<<iter<<std::endl;
+  for(int i=0;i<10;i++){
     result_o = Zero();
     t1=usecond();
     mCG(src_o,result_o);
@@ -170,24 +111,10 @@ int main (int argc, char ** argv)
     flops+= CGsiteflops*FrbGrid->gSites()*iters;
     std::cout << " SinglePrecision iterations/sec "<< iters/(t2-t1)*1000.*1000.<<std::endl;
     std::cout << " SinglePrecision GF/s "<< flops/(t2-t1)/1000.<<std::endl;
-    std::cout << " SinglePrecision error count "<< FlightRecorder::ErrorCount()<<std::endl;
-
-    assert(FlightRecorder::ErrorCount()==0);
-
-    std::cout << " FlightRecorder is OK! "<<std::endl;
-    iter ++;
-  } while (time(NULL) < (start + nsecs/10) );
-    
-  std::cout << GridLogMessage << "::::::::::::: Starting double precision CG" << std::endl;
+  }
+  std::cout << GridLogMessage << "::::::::::::: Starting regular CG" << std::endl;
   ConjugateGradient<LatticeFermionD> CG(1.0e-8,10000);
-  int i=0;
-  do { 
-    if ( i == 0 ) {
-      FlightRecorder::SetLoggingMode(FlightRecorder::LoggingModeRecord);
-    } else {
-      FlightRecorder::SetLoggingMode(FlightRecorder::LoggingModeVerify);
-    }
-    std::cerr << "******************* DOUBLE PRECISION SOLVE "<<i<<std::endl;
+  for(int i=0;i<1;i++){
     result_o_2 = Zero();
     t1=usecond();
     CG(HermOpEO,src_o,result_o_2);
@@ -195,20 +122,46 @@ int main (int argc, char ** argv)
     iters = CG.IterationsToComplete;
     flops = MdagMsiteflops*4*FrbGrid->gSites()*iters; 
     flops+= CGsiteflops*FrbGrid->gSites()*iters;
-
+    
     std::cout << " DoublePrecision iterations/sec "<< iters/(t2-t1)*1000.*1000.<<std::endl;
     std::cout << " DoublePrecision GF/s "<< flops/(t2-t1)/1000.<<std::endl;
-    std::cout << " DoublePrecision error count "<< FlightRecorder::ErrorCount()<<std::endl;
-    assert(FlightRecorder::ErrorCount()==0);
-    std::cout << " FlightRecorder is OK! "<<std::endl;
-    i++;
-  } while (time(NULL) < (start + nsecs) );
+  }
+  
+  //  MemoryManager::Print();
 
   LatticeFermionD diff_o(FrbGrid);
   RealD diff = axpy_norm(diff_o, -1.0, result_o, result_o_2);
 
   std::cout << GridLogMessage << "::::::::::::: Diff between mixed and regular CG: " << diff << std::endl;
-  assert(diff < 1e-4);
+
+  #ifdef HAVE_LIME
+  if( GridCmdOptionExists(argv,argv+argc,"--checksums") ){
   
+  std::string file1("./Propagator1");
+  emptyUserRecord record;
+  uint32_t nersc_csum;
+  uint32_t scidac_csuma;
+  uint32_t scidac_csumb;
+  typedef SpinColourVectorD   FermionD;
+  typedef vSpinColourVectorD vFermionD;
+
+  BinarySimpleMunger<FermionD,FermionD> munge;
+  std::string format = getFormatString<vFermionD>();
+  
+  BinaryIO::writeLatticeObject<vFermionD,FermionD>(result_o,file1,munge, 0, format,
+						   nersc_csum,scidac_csuma,scidac_csumb);
+
+  std::cout << GridLogMessage << " Mixed checksums "<<std::hex << scidac_csuma << " "<<scidac_csumb<<std::endl;
+
+  BinaryIO::writeLatticeObject<vFermionD,FermionD>(result_o_2,file1,munge, 0, format,
+						   nersc_csum,scidac_csuma,scidac_csumb);
+
+  std::cout << GridLogMessage << " CG checksums "<<std::hex << scidac_csuma << " "<<scidac_csumb<<std::endl;
+  }
+  #endif
+  }
+  
+  MemoryManager::Print();
+
   Grid_finalize();
 }