Debug code

.gitignore

@@ -1,7 +1,3 @@
-# Doxygen stuff
-html/*
-latex/*
-
 # Compiled Object files #
 #########################
 *.slo

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/Makefile.am

@@ -66,10 +66,6 @@ if BUILD_FERMION_REPS
   extra_sources+=$(ADJ_FERMION_FILES)
   extra_sources+=$(TWOIND_FERMION_FILES)
 endif
-if BUILD_SP
-    extra_sources+=$(SP_FERMION_FILES)
-    extra_sources+=$(SP_TWOIND_FERMION_FILES)
-endif
 
 lib_LIBRARIES = libGrid.a
 

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 -> A = (PRECISION) d -> A;
-  zd -> Delta = (ZOLO_PRECISION) d -> Delta;
+  zd -> Delta = (PRECISION) d -> Delta;
-  zd -> epsilon = (ZOLO_PRECISION) d -> epsilon;
+  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));
+  zd -> a = (PRECISION*) malloc(zd -> dn * sizeof(PRECISION));
-  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (ZOLO_PRECISION) d -> a[m];
+  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));
+  zd -> ap = (PRECISION*) malloc(zd -> dd * sizeof(PRECISION));
-  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (ZOLO_PRECISION) d -> ap[m];
+  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));
+  zd -> alpha = (PRECISION*) malloc(zd -> da * sizeof(PRECISION));
-  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (ZOLO_PRECISION) d -> alpha[m];
+  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));
+  zd -> beta = (PRECISION*) malloc(zd -> db * sizeof(PRECISION));
-  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (ZOLO_PRECISION) d -> beta[m];
+  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));
+  zd -> gamma = (PRECISION*) malloc(zd -> n * sizeof(PRECISION));
-  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (ZOLO_PRECISION) d -> gamma[m];
+  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 -> A = (PRECISION) d -> A;
-  zd -> Delta = (ZOLO_PRECISION) d -> Delta;
+  zd -> Delta = (PRECISION) d -> Delta;
-  zd -> epsilon = (ZOLO_PRECISION) d -> epsilon;
+  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));
+  zd -> a = (PRECISION*) malloc(zd -> dn * sizeof(PRECISION));
-  for (m = 0; m < zd -> dn; m++) zd -> a[m] = (ZOLO_PRECISION) d -> a[m];
+  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));
+  zd -> ap = (PRECISION*) malloc(zd -> dd * sizeof(PRECISION));
-  for (m = 0; m < zd -> dd; m++) zd -> ap[m] = (ZOLO_PRECISION) d -> ap[m];
+  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));
+  zd -> alpha = (PRECISION*) malloc(zd -> da * sizeof(PRECISION));
-  for (m = 0; m < zd -> da; m++) zd -> alpha[m] = (ZOLO_PRECISION) d -> alpha[m];
+  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));
+  zd -> beta = (PRECISION*) malloc(zd -> db * sizeof(PRECISION));
-  for (m = 0; m < zd -> db; m++) zd -> beta[m] = (ZOLO_PRECISION) d -> beta[m];
+  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));
+  zd -> gamma = (PRECISION*) malloc(zd -> n * sizeof(PRECISION));
-  for (m = 0; m < zd -> n; m++) zd -> gamma[m] = (ZOLO_PRECISION) d -> gamma[m];
+  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) 
+    ? higham((PRECISION) eps, n) 
-    : zolotarev((ZOLO_PRECISION) eps, n, type);
+    : 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
+#ifndef PRECISION
-#define ZOLO_PRECISION double
+#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* higham(PRECISION epsilon, int n) ;
-ZOLOTAREV_DATA* zolotarev(ZOLO_PRECISION epsilon, int n, int type);
+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/algorithms/iterative/ImplicitlyRestartedLanczos.h

@@ -419,15 +419,14 @@ until convergence
 	}
       }
 
-      if ( Nconv < Nstop ) {
+      if ( Nconv < Nstop )
 	std::cout << GridLogIRL << "Nconv ("<<Nconv<<") < Nstop ("<<Nstop<<")"<<std::endl;
-	std::cout << GridLogIRL << "returning Nstop vectors, the last "<< Nstop-Nconv << "of which might meet convergence criterion only approximately" <<std::endl;
+
-      }
       eval=eval2;
       
       //Keep only converged
-      eval.resize(Nstop);// was Nconv
+      eval.resize(Nconv);// Nstop?
-      evec.resize(Nstop,grid);// was Nconv
+      evec.resize(Nconv,grid);// Nstop?
       basisSortInPlace(evec,eval,reverse);
       
     }

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

@@ -27,7 +27,7 @@ Author: Christoph Lehner <christoph@lhnr.de>
 *************************************************************************************/
 /*  END LEGAL */
 
-#define Mheader "SharedMemoryMpi: "
+#define header "SharedMemoryMpi: "
 
 #include <Grid/GridCore.h>
 #include <pwd.h>
@@ -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>
@@ -176,8 +174,8 @@ void GlobalSharedMemory::Init(Grid_MPI_Comm comm)
   MPI_Comm_size(WorldShmComm     ,&WorldShmSize);
 
   if ( WorldRank == 0) {
-    std::cout << Mheader " World communicator of size " <<WorldSize << std::endl;  
+    std::cout << header " World communicator of size " <<WorldSize << std::endl;  
-    std::cout << Mheader " Node  communicator of size " <<WorldShmSize << std::endl;
+    std::cout << header " Node  communicator of size " <<WorldShmSize << std::endl;
   }
   // WorldShmComm, WorldShmSize, WorldShmRank
 
@@ -454,7 +452,7 @@ void GlobalSharedMemory::OptimalCommunicatorSharedMemory(const Coordinate &proce
 #ifdef GRID_MPI3_SHMGET
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
-  std::cout << Mheader "SharedMemoryAllocate "<< bytes<< " shmget implementation "<<std::endl;
+  std::cout << header "SharedMemoryAllocate "<< bytes<< " shmget implementation "<<std::endl;
   assert(_ShmSetup==1);
   assert(_ShmAlloc==0);
 
@@ -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 << header " 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,16 +574,13 @@ 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;
     exit(EXIT_FAILURE);  
   }
   if ( WorldRank == 0 ){
-    std::cout << WorldRank << Mheader " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
+    std::cout << WorldRank << header " SharedMemoryMPI.cc acceleratorAllocDevice "<< bytes 
 	      << "bytes at "<< std::hex<< ShmCommBuf << " - "<<(bytes-1+(uint64_t)ShmCommBuf) <<std::dec<<" for comms buffers " <<std::endl;
   }
   SharedMemoryZero(ShmCommBuf,bytes);
@@ -569,8 +604,8 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 #ifdef GRID_SYCL_LEVEL_ZERO_IPC
     typedef struct { int fd; pid_t pid ; ze_ipc_mem_handle_t ze; } clone_mem_t;
 
-    auto zeDevice    = cl::sycl::get_native<cl::sycl::backend::ext_oneapi_level_zero>(theGridAccelerator->get_device());
+    auto zeDevice    = cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_device());
-    auto zeContext   = cl::sycl::get_native<cl::sycl::backend::ext_oneapi_level_zero>(theGridAccelerator->get_context());
+    auto zeContext   = cl::sycl::get_native<cl::sycl::backend::level_zero>(theGridAccelerator->get_context());
       
     ze_ipc_mem_handle_t ihandle;
     clone_mem_t handle;
@@ -703,12 +738,13 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
   _ShmAllocBytes=bytes;
   _ShmAlloc=1;
 }
+#endif
 
 #else 
 #ifdef GRID_MPI3_SHMMMAP
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
-  std::cout << Mheader "SharedMemoryAllocate "<< bytes<< " MMAP implementation "<< GRID_SHM_PATH <<std::endl;
+  std::cout << header "SharedMemoryAllocate "<< bytes<< " MMAP implementation "<< GRID_SHM_PATH <<std::endl;
   assert(_ShmSetup==1);
   assert(_ShmAlloc==0);
   //////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -745,7 +781,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
     assert(((uint64_t)ptr&0x3F)==0);
     close(fd);
     WorldShmCommBufs[r] =ptr;
-    //    std::cout << Mheader "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
+    //    std::cout << header "Set WorldShmCommBufs["<<r<<"]="<<ptr<< "("<< bytes<< "bytes)"<<std::endl;
   }
   _ShmAlloc=1;
   _ShmAllocBytes  = bytes;
@@ -755,7 +791,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 #ifdef GRID_MPI3_SHM_NONE
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 {
-  std::cout << Mheader "SharedMemoryAllocate "<< bytes<< " MMAP anonymous implementation "<<std::endl;
+  std::cout << header "SharedMemoryAllocate "<< bytes<< " MMAP anonymous implementation "<<std::endl;
   assert(_ShmSetup==1);
   assert(_ShmAlloc==0);
   //////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -802,7 +838,7 @@ void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 ////////////////////////////////////////////////////////////////////////////////////////////
 void GlobalSharedMemory::SharedMemoryAllocate(uint64_t bytes, int flags)
 { 
-  std::cout << Mheader "SharedMemoryAllocate "<< bytes<< " SHMOPEN implementation "<<std::endl;
+  std::cout << header "SharedMemoryAllocate "<< bytes<< " SHMOPEN implementation "<<std::endl;
   assert(_ShmSetup==1);
   assert(_ShmAlloc==0); 
   MPI_Barrier(WorldShmComm);
@@ -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/cshift/Cshift_common.h

@@ -29,27 +29,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 
 NAMESPACE_BEGIN(Grid);
 
-extern std::vector<std::pair<int,int> > Cshift_table; 
+extern Vector<std::pair<int,int> > Cshift_table; 
-extern commVector<std::pair<int,int> > Cshift_table_device; 
 
-inline std::pair<int,int> *MapCshiftTable(void)
-{
-  // GPU version
-#ifdef ACCELERATOR_CSHIFT    
-  uint64_t sz=Cshift_table.size();
-  if (Cshift_table_device.size()!=sz )    {
-    Cshift_table_device.resize(sz);
-  }
-  acceleratorCopyToDevice((void *)&Cshift_table[0],
-			  (void *)&Cshift_table_device[0],
-			  sizeof(Cshift_table[0])*sz);
-
-  return &Cshift_table_device[0];
-#else 
-  return &Cshift_table[0];
-#endif
-  // CPU version use identify map
-}
 ///////////////////////////////////////////////////////////////////
 // Gather for when there is no need to SIMD split 
 ///////////////////////////////////////////////////////////////////
@@ -93,8 +74,8 @@ Gather_plane_simple (const Lattice<vobj> &rhs,cshiftVector<vobj> &buffer,int dim
   }
   {
     auto buffer_p = & buffer[0];
-    auto table = MapCshiftTable();
+    auto table = &Cshift_table[0];
-#ifdef ACCELERATOR_CSHIFT
+#ifdef ACCELERATOR_CSHIFT    
     autoView(rhs_v , rhs, AcceleratorRead);
     accelerator_for(i,ent,vobj::Nsimd(),{
 	coalescedWrite(buffer_p[table[i].first],coalescedRead(rhs_v[table[i].second]));
@@ -244,7 +225,7 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,cshiftVector<
   
   {
     auto buffer_p = & buffer[0];
-    auto table = MapCshiftTable();
+    auto table = &Cshift_table[0];
 #ifdef ACCELERATOR_CSHIFT    
     autoView( rhs_v, rhs, AcceleratorWrite);
     accelerator_for(i,ent,vobj::Nsimd(),{
@@ -316,6 +297,30 @@ template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerA
   }
 }
 
+#if (defined(GRID_CUDA) || defined(GRID_HIP)) && defined(ACCELERATOR_CSHIFT)
+
+template <typename T>
+T iDivUp(T a, T b) // Round a / b to nearest higher integer value
+{ return (a % b != 0) ? (a / b + 1) : (a / b); }
+
+template <typename T>
+__global__ void populate_Cshift_table(T* vector, T lo, T ro, T e1, T e2, T stride)
+{
+    int idx = blockIdx.x*blockDim.x + threadIdx.x;
+    if (idx >= e1*e2) return;
+
+    int n, b, o;
+
+    n = idx / e2;
+    b = idx % e2;
+    o = n*stride + b;
+
+    vector[2*idx + 0] = lo + o;
+    vector[2*idx + 1] = ro + o;
+}
+
+#endif
+
 //////////////////////////////////////////////////////
 // local to node block strided copies
 //////////////////////////////////////////////////////
@@ -340,12 +345,20 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs
   int ent=0;
 
   if(cbmask == 0x3 ){
+#if (defined(GRID_CUDA) || defined(GRID_HIP)) && defined(ACCELERATOR_CSHIFT)
+    ent = e1*e2;
+    dim3 blockSize(acceleratorThreads());
+    dim3 gridSize(iDivUp((unsigned int)ent, blockSize.x));
+    populate_Cshift_table<<<gridSize, blockSize>>>(&Cshift_table[0].first, lo, ro, e1, e2, stride);
+    accelerator_barrier();
+#else
     for(int n=0;n<e1;n++){
       for(int b=0;b<e2;b++){
         int o =n*stride+b;
 	Cshift_table[ent++] = std::pair<int,int>(lo+o,ro+o);
       }
     }
+#endif
   } else { 
     for(int n=0;n<e1;n++){
       for(int b=0;b<e2;b++){
@@ -359,7 +372,7 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs
   }
 
   {
-    auto table = MapCshiftTable();
+    auto table = &Cshift_table[0];
 #ifdef ACCELERATOR_CSHIFT    
     autoView(rhs_v , rhs, AcceleratorRead);
     autoView(lhs_v , lhs, AcceleratorWrite);
@@ -396,11 +409,19 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vo
   int ent=0;
 
   if ( cbmask == 0x3 ) {
+#if (defined(GRID_CUDA) || defined(GRID_HIP)) && defined(ACCELERATOR_CSHIFT)
+    ent = e1*e2;
+    dim3 blockSize(acceleratorThreads());
+    dim3 gridSize(iDivUp((unsigned int)ent, blockSize.x));
+    populate_Cshift_table<<<gridSize, blockSize>>>(&Cshift_table[0].first, lo, ro, e1, e2, stride);
+    accelerator_barrier();
+#else
     for(int n=0;n<e1;n++){
     for(int b=0;b<e2;b++){
       int o  =n*stride;
       Cshift_table[ent++] = std::pair<int,int>(lo+o+b,ro+o+b);
     }}
+#endif
   } else {
     for(int n=0;n<e1;n++){
     for(int b=0;b<e2;b++){
@@ -411,7 +432,7 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vo
   }
 
   {
-    auto table = MapCshiftTable();
+    auto table = &Cshift_table[0];
 #ifdef ACCELERATOR_CSHIFT    
     autoView( rhs_v, rhs, AcceleratorRead);
     autoView( lhs_v, lhs, AcceleratorWrite);

Grid/cshift/Cshift_mpi.h

@@ -52,8 +52,7 @@ template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension
   int comm_dim        = rhs.Grid()->_processors[dimension] >1 ;
   int splice_dim      = rhs.Grid()->_simd_layout[dimension]>1 && (comm_dim);
 
-  RealD t1,t0;
+
-  t0=usecond();
   if ( !comm_dim ) {
     //std::cout << "CSHIFT: Cshift_local" <<std::endl;
     Cshift_local(ret,rhs,dimension,shift); // Handles checkerboarding
@@ -64,8 +63,6 @@ template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension
     //std::cout << "CSHIFT: Cshift_comms" <<std::endl;
     Cshift_comms(ret,rhs,dimension,shift);
   }
-  t1=usecond();
-  //  std::cout << GridLogPerformance << "Cshift took "<< (t1-t0)/1e3 << " ms"<<std::endl;
   return ret;
 }
 
@@ -130,20 +127,16 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
     
   int cb= (cbmask==0x2)? Odd : Even;
   int sshift= rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,cb);
-  RealD tcopy=0.0;
+
-  RealD tgather=0.0;
-  RealD tscatter=0.0;
-  RealD tcomms=0.0;
-  uint64_t xbytes=0;
   for(int x=0;x<rd;x++){       
 
     int sx        =  (x+sshift)%rd;
     int comm_proc = ((x+sshift)/rd)%pd;
     
     if (comm_proc==0) {
-      tcopy-=usecond();
+
       Copy_plane(ret,rhs,dimension,x,sx,cbmask); 
-      tcopy+=usecond();
+
     } else {
 
       int words = buffer_size;
@@ -151,39 +144,26 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
 
       int bytes = words * sizeof(vobj);
 
-      tgather-=usecond();
       Gather_plane_simple (rhs,send_buf,dimension,sx,cbmask);
-      tgather+=usecond();
 
       //      int rank           = grid->_processor;
       int recv_from_rank;
       int xmit_to_rank;
       grid->ShiftedRanks(dimension,comm_proc,xmit_to_rank,recv_from_rank);
-      
+
-      tcomms-=usecond();
+      grid->Barrier();
-      //      grid->Barrier();
 
       grid->SendToRecvFrom((void *)&send_buf[0],
 			   xmit_to_rank,
 			   (void *)&recv_buf[0],
 			   recv_from_rank,
 			   bytes);
-      xbytes+=bytes;
-      //      grid->Barrier();
-      tcomms+=usecond();
 
-      tscatter-=usecond();
+      grid->Barrier();
+
       Scatter_plane_simple (ret,recv_buf,dimension,x,cbmask);
-      tscatter+=usecond();
     }
   }
-  /*
-  std::cout << GridLogPerformance << " Cshift copy    "<<tcopy/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift gather  "<<tgather/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift scatter "<<tscatter/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift comm    "<<tcomms/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift BW      "<<(2.0*xbytes)/tcomms<<" MB/s "<<2*xbytes<< " Bytes "<<std::endl;
-  */
 }
 
 template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
@@ -210,12 +190,6 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
   assert(shift>=0);
   assert(shift<fd);
 
-  RealD tcopy=0.0;
-  RealD tgather=0.0;
-  RealD tscatter=0.0;
-  RealD tcomms=0.0;
-  uint64_t xbytes=0;
-  
   int permute_type=grid->PermuteType(dimension);
 
   ///////////////////////////////////////////////
@@ -253,9 +227,7 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
       pointers[i] = &send_buf_extract[i][0];
     }
     int sx   = (x+sshift)%rd;
-    tgather-=usecond();
     Gather_plane_extract(rhs,pointers,dimension,sx,cbmask);
-    tgather+=usecond();
 
     for(int i=0;i<Nsimd;i++){
       
@@ -280,8 +252,7 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
       if(nbr_proc){
 	grid->ShiftedRanks(dimension,nbr_proc,xmit_to_rank,recv_from_rank); 
 
-	tcomms-=usecond();
+	grid->Barrier();
-	//	grid->Barrier();
 
 	send_buf_extract_mpi = &send_buf_extract[nbr_lane][0];
 	recv_buf_extract_mpi = &recv_buf_extract[i][0];
@@ -291,9 +262,7 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
 			     recv_from_rank,
 			     bytes);
 
-	xbytes+=bytes;
+	grid->Barrier();
-	//	grid->Barrier();
-	tcomms+=usecond();
 
 	rpointers[i] = &recv_buf_extract[i][0];
       } else { 
@@ -301,17 +270,9 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
       }
 
     }
-    tscatter-=usecond();
     Scatter_plane_merge(ret,rpointers,dimension,x,cbmask);
-    tscatter+=usecond();
   }
-  /*
+
-  std::cout << GridLogPerformance << " Cshift (s) copy    "<<tcopy/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) gather  "<<tgather/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) scatter "<<tscatter/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) comm    "<<tcomms/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift BW      "<<(2.0*xbytes)/tcomms<<" MB/s "<<2*xbytes<< " Bytes "<<std::endl;
-  */
 }
 #else 
 template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
@@ -331,11 +292,6 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
   assert(comm_dim==1);
   assert(shift>=0);
   assert(shift<fd);
-  RealD tcopy=0.0;
-  RealD tgather=0.0;
-  RealD tscatter=0.0;
-  RealD tcomms=0.0;
-  uint64_t xbytes=0;
   
   int buffer_size = rhs.Grid()->_slice_nblock[dimension]*rhs.Grid()->_slice_block[dimension];
   static cshiftVector<vobj> send_buf_v; send_buf_v.resize(buffer_size);
@@ -359,9 +315,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
     
     if (comm_proc==0) {
 
-      tcopy-=usecond();
       Copy_plane(ret,rhs,dimension,x,sx,cbmask); 
-      tcopy+=usecond();
 
     } else {
 
@@ -370,9 +324,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
 
       int bytes = words * sizeof(vobj);
 
-      tgather-=usecond();
       Gather_plane_simple (rhs,send_buf_v,dimension,sx,cbmask);
-      tgather+=usecond();
 
       //      int rank           = grid->_processor;
       int recv_from_rank;
@@ -380,8 +332,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
       grid->ShiftedRanks(dimension,comm_proc,xmit_to_rank,recv_from_rank);
 
 
-      tcomms-=usecond();
+      grid->Barrier();
-      //      grid->Barrier();
 
       acceleratorCopyDeviceToDevice((void *)&send_buf_v[0],(void *)&send_buf[0],bytes);
       grid->SendToRecvFrom((void *)&send_buf[0],
@@ -389,24 +340,13 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
 			   (void *)&recv_buf[0],
 			   recv_from_rank,
 			   bytes);
-      xbytes+=bytes;
       acceleratorCopyDeviceToDevice((void *)&recv_buf[0],(void *)&recv_buf_v[0],bytes);
 
-      //      grid->Barrier();
+      grid->Barrier();
-      tcomms+=usecond();
 
-      tscatter-=usecond();
       Scatter_plane_simple (ret,recv_buf_v,dimension,x,cbmask);
-      tscatter+=usecond();
     }
   }
-  /*
-  std::cout << GridLogPerformance << " Cshift copy    "<<tcopy/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift gather  "<<tgather/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift scatter "<<tscatter/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift comm    "<<tcomms/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift BW      "<<(2.0*xbytes)/tcomms<<" MB/s "<<2*xbytes<< " Bytes "<<std::endl;
-  */
 }
 
 template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
@@ -432,11 +372,6 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
   assert(simd_layout==2);
   assert(shift>=0);
   assert(shift<fd);
-  RealD tcopy=0.0;
-  RealD tgather=0.0;
-  RealD tscatter=0.0;
-  RealD tcomms=0.0;
-  uint64_t xbytes=0;
 
   int permute_type=grid->PermuteType(dimension);
 
@@ -479,10 +414,8 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
     for(int i=0;i<Nsimd;i++){       
       pointers[i] = &send_buf_extract[i][0];
     }
-    tgather-=usecond();
     int sx   = (x+sshift)%rd;
     Gather_plane_extract(rhs,pointers,dimension,sx,cbmask);
-    tgather+=usecond();
 
     for(int i=0;i<Nsimd;i++){
       
@@ -507,8 +440,7 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
       if(nbr_proc){
 	grid->ShiftedRanks(dimension,nbr_proc,xmit_to_rank,recv_from_rank); 
 
-	tcomms-=usecond();
+	grid->Barrier();
-	//	grid->Barrier();
 
 	acceleratorCopyDeviceToDevice((void *)&send_buf_extract[nbr_lane][0],(void *)send_buf_extract_mpi,bytes);
 	grid->SendToRecvFrom((void *)send_buf_extract_mpi,
@@ -517,28 +449,17 @@ template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vo
 			     recv_from_rank,
 			     bytes);
 	acceleratorCopyDeviceToDevice((void *)recv_buf_extract_mpi,(void *)&recv_buf_extract[i][0],bytes);
-	xbytes+=bytes;
 
-	//	grid->Barrier();
+	grid->Barrier();
-	tcomms+=usecond();
 	rpointers[i] = &recv_buf_extract[i][0];
       } else { 
 	rpointers[i] = &send_buf_extract[nbr_lane][0];
       }
 
     }
-    tscatter-=usecond();
     Scatter_plane_merge(ret,rpointers,dimension,x,cbmask);
-    tscatter+=usecond();
-
   }
-  /*
+
-  std::cout << GridLogPerformance << " Cshift (s) copy    "<<tcopy/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) gather  "<<tgather/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) scatter "<<tscatter/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift (s) comm    "<<tcomms/1e3<<" ms"<<std::endl;
-  std::cout << GridLogPerformance << " Cshift BW      "<<(2.0*xbytes)/tcomms<<" MB/s"<<std::endl;
-  */
 }
 #endif
 NAMESPACE_END(Grid); 

Grid/cshift/Cshift_table.cc

@@ -1,5 +1,4 @@
 #include <Grid/GridCore.h>       
 NAMESPACE_BEGIN(Grid);
-std::vector<std::pair<int,int> > Cshift_table; 
+Vector<std::pair<int,int> > Cshift_table; 
-commVector<std::pair<int,int> > Cshift_table_device; 
 NAMESPACE_END(Grid);

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,4 @@ 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/PaddedCell.h>
+#include <Grid/lattice/Lattice_crc.h>

Grid/lattice/Lattice_ET.h

@@ -345,9 +345,7 @@ GridUnopClass(UnaryNot, Not(a));
 GridUnopClass(UnaryTrace, trace(a));
 GridUnopClass(UnaryTranspose, transpose(a));
 GridUnopClass(UnaryTa, Ta(a));
-GridUnopClass(UnarySpTa, SpTa(a));
 GridUnopClass(UnaryProjectOnGroup, ProjectOnGroup(a));
-GridUnopClass(UnaryProjectOnSpGroup, ProjectOnSpGroup(a));
 GridUnopClass(UnaryTimesI, timesI(a));
 GridUnopClass(UnaryTimesMinusI, timesMinusI(a));
 GridUnopClass(UnaryAbs, abs(a));
@@ -458,9 +456,7 @@ GRID_DEF_UNOP(operator!, UnaryNot);
 GRID_DEF_UNOP(trace, UnaryTrace);
 GRID_DEF_UNOP(transpose, UnaryTranspose);
 GRID_DEF_UNOP(Ta, UnaryTa);
-GRID_DEF_UNOP(SpTa, UnarySpTa);
 GRID_DEF_UNOP(ProjectOnGroup, UnaryProjectOnGroup);
-GRID_DEF_UNOP(ProjectOnSpGroup, UnaryProjectOnSpGroup);
 GRID_DEF_UNOP(timesI, UnaryTimesI);
 GRID_DEF_UNOP(timesMinusI, UnaryTimesMinusI);
 GRID_DEF_UNOP(abs, UnaryAbs);  // abs overloaded in cmath C++98; DON'T do the

Grid/lattice/Lattice_arith.h

@@ -270,42 +270,5 @@ RealD axpby_norm(Lattice<vobj> &ret,sobj a,sobj b,const Lattice<vobj> &x,const L
     return axpby_norm_fast(ret,a,b,x,y);
 }
 
-/// Trace product
-template<class obj> auto traceProduct(const Lattice<obj> &rhs_1,const Lattice<obj> &rhs_2)
-  -> Lattice<decltype(trace(obj()))>
-{
-  typedef decltype(trace(obj())) robj;
-  Lattice<robj> ret_i(rhs_1.Grid());
-  autoView( rhs1 , rhs_1, AcceleratorRead);
-  autoView( rhs2 , rhs_2, AcceleratorRead);
-  autoView( ret , ret_i, AcceleratorWrite);
-  ret.Checkerboard() = rhs_1.Checkerboard();
-  accelerator_for(ss,rhs1.size(),obj::Nsimd(),{
-      coalescedWrite(ret[ss],traceProduct(rhs1(ss),rhs2(ss)));
-  });
-  return ret_i;
-}
-
-template<class obj1,class obj2> auto traceProduct(const Lattice<obj1> &rhs_1,const obj2 &rhs2)
-  -> Lattice<decltype(trace(obj1()))>
-{
-  typedef decltype(trace(obj1())) robj;
-  Lattice<robj> ret_i(rhs_1.Grid());
-  autoView( rhs1 , rhs_1, AcceleratorRead);
-  autoView( ret , ret_i, AcceleratorWrite);
-  ret.Checkerboard() = rhs_1.Checkerboard();
-  accelerator_for(ss,rhs1.size(),obj1::Nsimd(),{
-      coalescedWrite(ret[ss],traceProduct(rhs1(ss),rhs2));
-  });
-  return ret_i;
-}
-template<class obj1,class obj2> auto traceProduct(const obj2 &rhs_2,const Lattice<obj1> &rhs_1)
-  -> Lattice<decltype(trace(obj1()))>
-{
-  return traceProduct(rhs_1,rhs_2);
-}
-
-
-
 NAMESPACE_END(Grid);
 #endif

Grid/lattice/Lattice_basis.h

@@ -62,7 +62,7 @@ void basisRotate(VField &basis,Matrix& Qt,int j0, int j1, int k0,int k1,int Nm)
     basis_v.push_back(basis[k].View(AcceleratorWrite));
   }
 
-#if ( !(defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)) )
+#if ( (!defined(GRID_CUDA)) )
   int max_threads = thread_max();
   Vector < vobj > Bt(Nm * max_threads);
   thread_region

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

@@ -31,7 +31,6 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #if defined(GRID_SYCL)
 #include <Grid/lattice/Lattice_reduction_sycl.h>
 #endif
-#include <Grid/lattice/Lattice_slicesum_core.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -281,29 +280,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;
 }
 
@@ -467,10 +448,19 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   int e1=    grid->_slice_nblock[orthogdim];
   int e2=    grid->_slice_block [orthogdim];
   int stride=grid->_slice_stride[orthogdim];
-  int ostride=grid->_ostride[orthogdim];
+
-  
+  // sum over reduced dimension planes, breaking out orthog dir
-  //Reduce Data down to lvSum
+  // Parallel over orthog direction
-  sliceSumReduction(Data,lvSum,rd, e1,e2,stride,ostride,Nsimd);
+  autoView( Data_v, Data, CpuRead);
+  thread_for( r,rd, {
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+	int ss= so+n*stride+b;
+	lvSum[r]=lvSum[r]+Data_v[ss];
+      }
+    }
+  });
 
   // Sum across simd lanes in the plane, breaking out orthog dir.
   Coordinate icoor(Nd);
@@ -514,7 +504,6 @@ sliceSum(const Lattice<vobj> &Data,int orthogdim)
   return result;
 }
 
-
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {

Grid/lattice/Lattice_reduction_gpu.h

@@ -30,7 +30,7 @@ int getNumBlocksAndThreads(const Iterator n, const size_t sizeofsobj, Iterator &
   cudaGetDevice(&device);
 #endif
 #ifdef GRID_HIP
-  auto r=hipGetDevice(&device);
+  hipGetDevice(&device);
 #endif
   
   Iterator warpSize            = gpu_props[device].warpSize;

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

@@ -152,7 +152,6 @@ public:
 #ifdef RNG_FAST_DISCARD
   static void Skip(RngEngine &eng,uint64_t site)
   {
-#if 0
     /////////////////////////////////////////////////////////////////////////////////////
     // Skip by 2^40 elements between successive lattice sites
     // This goes by 10^12.
@@ -163,9 +162,9 @@ public:
     // tens of seconds per trajectory so this is clean in all reasonable cases,
     // and margin of safety is orders of magnitude.
     // We could hack Sitmo to skip in the higher order words of state if necessary
-    //
+      //
-    // Replace with 2^30 ; avoid problem on large volumes
+      // Replace with 2^30 ; avoid problem on large volumes
-    //
+      //
     /////////////////////////////////////////////////////////////////////////////////////
     //      uint64_t skip = site+1;  //   Old init Skipped then drew.  Checked compat with faster init
     const int shift = 30;
@@ -180,9 +179,6 @@ public:
     assert((skip >> shift)==site); // check for overflow
 
     eng.discard(skip);
-#else
-    eng.discardhi(site);
-#endif
     //      std::cout << " Engine  " <<site << " state " <<eng<<std::endl;
   } 
 #endif
@@ -411,7 +407,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 +424,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 +445,29 @@ public:
 	
 	int l_idx=generator_idx(o_idx,i_idx);
 	_generators[l_idx] = master_engine;
-	if ( britney ) { 
+	Skip(_generators[l_idx],gidx); // Skip to next RNG sequence
-	  Skip(_generators[l_idx],l_idx); // Skip to next RNG sequence
+    });
-	} else { 	
+#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/lattice/Lattice_slicesum_core.h

@@ -1,224 +0,0 @@
-#pragma once
-
-#if defined(GRID_CUDA)
-
-#include <cub/cub.cuh>
-#define gpucub cub
-#define gpuError_t cudaError_t
-#define gpuSuccess cudaSuccess
-
-#elif defined(GRID_HIP)
-
-#include <hipcub/hipcub.hpp>
-#define gpucub hipcub
-#define gpuError_t hipError_t
-#define gpuSuccess hipSuccess
-
-#endif
-
-
-NAMESPACE_BEGIN(Grid);
-
-
-#if defined(GRID_CUDA) || defined(GRID_HIP)
-template<class vobj> inline void sliceSumReduction_cub_small(const vobj *Data, Vector<vobj> &lvSum, const int rd, const int e1, const int e2, const int stride, const int ostride, const int Nsimd) {
-  size_t subvol_size = e1*e2;
-  commVector<vobj> reduction_buffer(rd*subvol_size);
-  auto rb_p = &reduction_buffer[0];
-  vobj zero_init;
-  zeroit(zero_init);
-
-  
-  void *temp_storage_array = NULL;
-  size_t temp_storage_bytes = 0;
-  vobj *d_out;
-  int* d_offsets;
-
-  std::vector<int> offsets(rd+1,0);
-
-  for (int i = 0; i < offsets.size(); i++) {
-    offsets[i] = i*subvol_size;
-  }
-  
-  //Allocate memory for output and offset arrays on device
-  d_out = static_cast<vobj*>(acceleratorAllocDevice(rd*sizeof(vobj)));
-  
-  d_offsets = static_cast<int*>(acceleratorAllocDevice((rd+1)*sizeof(int)));
-  
-  //copy offsets to device
-  acceleratorCopyToDeviceAsync(&offsets[0],d_offsets,sizeof(int)*(rd+1),computeStream);
-  
-  
-  gpuError_t gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p,d_out, rd, d_offsets, d_offsets+1, ::gpucub::Sum(), zero_init, computeStream);
-  if (gpuErr!=gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce (setup)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  //allocate memory for temp_storage_array  
-  temp_storage_array = acceleratorAllocDevice(temp_storage_bytes);
-  
-  //prepare buffer for reduction
-  //use non-blocking accelerator_for to avoid syncs (ok because we submit to same computeStream)
-  //use 2d accelerator_for to avoid launch latencies found when serially looping over rd 
-  accelerator_for2dNB( s,subvol_size, r,rd, Nsimd,{ 
-  
-    int n = s / e2;
-    int b = s % e2;
-    int so=r*ostride; // base offset for start of plane 
-    int ss= so+n*stride+b;
-
-    coalescedWrite(rb_p[r*subvol_size+s], coalescedRead(Data[ss]));
-
-  });
-  
-  //issue segmented reductions in computeStream
-  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p, d_out, rd, d_offsets, d_offsets+1,::gpucub::Sum(), zero_init, computeStream);
-  if (gpuErr!=gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-  
-  acceleratorCopyFromDeviceAsync(d_out,&lvSum[0],rd*sizeof(vobj),computeStream);
-  
-  //sync after copy
-  accelerator_barrier();
- 
-  acceleratorFreeDevice(temp_storage_array);
-  acceleratorFreeDevice(d_out);
-  acceleratorFreeDevice(d_offsets);
-  
-
-}
-#endif 
-
-
-#if defined(GRID_SYCL)
-template<class vobj> inline void sliceSumReduction_sycl_small(const vobj *Data, Vector <vobj> &lvSum, const int  &rd, const int &e1, const int &e2, const int &stride, const int &ostride, const int &Nsimd)
-{
-  size_t subvol_size = e1*e2;
-
-  vobj *mysum = (vobj *) malloc_shared(rd*sizeof(vobj),*theGridAccelerator);
-  vobj vobj_zero;
-  zeroit(vobj_zero);
-  for (int r = 0; r<rd; r++) { 
-    mysum[r] = vobj_zero; 
-  }
-
-  commVector<vobj> reduction_buffer(rd*subvol_size);    
-
-  auto rb_p = &reduction_buffer[0];
-
-  // autoView(Data_v, Data, AcceleratorRead);
-
-  //prepare reduction buffer 
-  accelerator_for2d( s,subvol_size, r,rd, (size_t)Nsimd,{ 
-  
-      int n = s / e2;
-      int b = s % e2;
-      int so=r*ostride; // base offset for start of plane 
-      int ss= so+n*stride+b;
-
-      coalescedWrite(rb_p[r*subvol_size+s], coalescedRead(Data[ss]));
-
-  });
-
-  for (int r = 0; r < rd; r++) {
-      theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
-          auto Reduction = cl::sycl::reduction(&mysum[r],std::plus<>());
-          cgh.parallel_for(cl::sycl::range<1>{subvol_size},
-          Reduction,
-          [=](cl::sycl::id<1> item, auto &sum) {
-              auto s = item[0];
-              sum += rb_p[r*subvol_size+s];
-          });
-      });
-      
-     
-  }
-  theGridAccelerator->wait();
-  for (int r = 0; r < rd; r++) {
-    lvSum[r] = mysum[r];
-  }
-  free(mysum,*theGridAccelerator);
-}
-#endif
-
-template<class vobj> inline void sliceSumReduction_large(const vobj *Data, Vector<vobj> &lvSum, const int rd, const int e1, const int e2, const int stride, const int ostride, const int Nsimd) {
-  typedef typename vobj::vector_type vector;
-  const int words = sizeof(vobj)/sizeof(vector);
-  const int osites = rd*e1*e2;
-  commVector<vector>buffer(osites);
-  vector *dat = (vector *)Data;
-  vector *buf = &buffer[0];
-  Vector<vector> lvSum_small(rd);
-  vector *lvSum_ptr = (vector *)&lvSum[0];
-
-  for (int w = 0; w < words; w++) {
-    accelerator_for(ss,osites,1,{
-	    buf[ss] = dat[ss*words+w];
-    });
-
-    #if defined(GRID_CUDA) || defined(GRID_HIP)
-      sliceSumReduction_cub_small(buf,lvSum_small,rd,e1,e2,stride, ostride,Nsimd);
-    #elif defined(GRID_SYCL)
-      sliceSumReduction_sycl_small(buf,lvSum_small,rd,e1,e2,stride, ostride,Nsimd);
-    #endif
-
-    for (int r = 0; r < rd; r++) {
-      lvSum_ptr[w+words*r]=lvSum_small[r];
-    }
-
-  }
-
-  
-}
-
-template<class vobj> inline void sliceSumReduction_gpu(const Lattice<vobj> &Data, Vector<vobj> &lvSum, const int rd, const int e1, const int e2, const int stride, const int ostride, const int Nsimd)
-{
-  autoView(Data_v, Data, AcceleratorRead); //reduction libraries cannot deal with large vobjs so we split into small/large case.
-    if constexpr (sizeof(vobj) <= 256) { 
-
-      #if defined(GRID_CUDA) || defined(GRID_HIP)
-        sliceSumReduction_cub_small(&Data_v[0], lvSum, rd, e1, e2, stride, ostride, Nsimd);
-      #elif defined (GRID_SYCL)
-        sliceSumReduction_sycl_small(&Data_v[0], lvSum, rd, e1, e2, stride, ostride, Nsimd);
-      #endif
-
-    }
-    else {
-      sliceSumReduction_large(&Data_v[0], lvSum, rd, e1, e2, stride, ostride, Nsimd);
-    }
-}
-
-
-template<class vobj> inline void sliceSumReduction_cpu(const Lattice<vobj> &Data, Vector<vobj> &lvSum, const int &rd, const int &e1, const int &e2, const int &stride, const int &ostride, const int &Nsimd)
-{
-  // sum over reduced dimension planes, breaking out orthog dir
-  // Parallel over orthog direction
-  autoView( Data_v, Data, CpuRead);
-  thread_for( r,rd, {
-    int so=r*ostride; // base offset for start of plane 
-    for(int n=0;n<e1;n++){
-      for(int b=0;b<e2;b++){
-        int ss= so+n*stride+b;
-        lvSum[r]=lvSum[r]+Data_v[ss];
-      }
-    }
-  });
-}
-
-template<class vobj> inline void sliceSumReduction(const Lattice<vobj> &Data, Vector<vobj> &lvSum, const int &rd, const int &e1, const int &e2, const int &stride, const int &ostride, const int &Nsimd) 
-{
-  #if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
-  
-  sliceSumReduction_gpu(Data, lvSum, rd, e1, e2, stride, ostride, Nsimd);
-  
-  #else
-  sliceSumReduction_cpu(Data, lvSum, rd, e1, e2, stride, ostride, Nsimd);
-
-  #endif
-}
-
-
-NAMESPACE_END(Grid);

Grid/lattice/Lattice_trace.h

@@ -66,65 +66,6 @@ inline auto TraceIndex(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<
   return ret;
 };
 
-template<int N, class Vec>
-Lattice<iScalar<iScalar<iScalar<Vec> > > > Determinant(const Lattice<iScalar<iScalar<iMatrix<Vec, N> > > > &Umu)
-{
-  GridBase *grid=Umu.Grid();
-  auto lvol = grid->lSites();
-  Lattice<iScalar<iScalar<iScalar<Vec> > > > ret(grid);
-  typedef typename Vec::scalar_type scalar;
-  autoView(Umu_v,Umu,CpuRead);
-  autoView(ret_v,ret,CpuWrite);
-  thread_for(site,lvol,{
-    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
-    Coordinate lcoor;
-    grid->LocalIndexToLocalCoor(site, lcoor);
-    iScalar<iScalar<iMatrix<scalar, N> > > Us;
-    peekLocalSite(Us, Umu_v, lcoor);
-    for(int i=0;i<N;i++){
-      for(int j=0;j<N;j++){
-	scalar tmp= Us()()(i,j);
-	ComplexD ztmp(real(tmp),imag(tmp));
-	EigenU(i,j)=ztmp;
-      }}
-    ComplexD detD  = EigenU.determinant();
-    typename Vec::scalar_type det(detD.real(),detD.imag());
-    pokeLocalSite(det,ret_v,lcoor);
-  });
-  return ret;
-}
-
-template<int N>
-Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
-{
-  GridBase *grid=Umu.Grid();
-  auto lvol = grid->lSites();
-  Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
-  
-  autoView(Umu_v,Umu,CpuRead);
-  autoView(ret_v,ret,CpuWrite);
-  thread_for(site,lvol,{
-    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
-    Coordinate lcoor;
-    grid->LocalIndexToLocalCoor(site, lcoor);
-    iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
-    iScalar<iScalar<iMatrix<ComplexD, N> > > Ui;
-    peekLocalSite(Us, Umu_v, lcoor);
-    for(int i=0;i<N;i++){
-      for(int j=0;j<N;j++){
-	EigenU(i,j) = Us()()(i,j);
-      }}
-    Eigen::MatrixXcd EigenUinv = EigenU.inverse();
-    for(int i=0;i<N;i++){
-      for(int j=0;j<N;j++){
-	Ui()()(i,j) = EigenUinv(i,j);
-      }}
-    pokeLocalSite(Ui,ret_v,lcoor);
-  });
-  return ret;
-}
-
-
 NAMESPACE_END(Grid);
 #endif
 

Grid/lattice/Lattice_transfer.h

@@ -469,13 +469,15 @@ inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
   Coordinate fine_rdimensions = fine->_rdimensions;
   Coordinate coarse_rdimensions = coarse->_rdimensions;
 
+  vobj zz = Zero();
+  
   accelerator_for(sc,coarse->oSites(),1,{
 
       // One thread per sub block
       Coordinate coor_c(_ndimension);
       Lexicographic::CoorFromIndex(coor_c,sc,coarse_rdimensions);  // Block coordinate
 
-      vobj cd = Zero();
+      vobj cd = zz;
       
       for(int sb=0;sb<blockVol;sb++){
 
@@ -695,68 +697,8 @@ void localCopyRegion(const Lattice<vobj> &From,Lattice<vobj> & To,Coordinate Fro
   for(int d=0;d<nd;d++){
     assert(Fg->_processors[d]  == Tg->_processors[d]);
   }
+
   // the above should guarantee that the operations are local
-  
-#if 1
-
-  size_t nsite = 1;
-  for(int i=0;i<nd;i++) nsite *= RegionSize[i];
-  
-  size_t tbytes = 4*nsite*sizeof(int);
-  int *table = (int*)malloc(tbytes);
- 
-  thread_for(idx, nsite, {
-      Coordinate from_coor, to_coor;
-      size_t rem = idx;
-      for(int i=0;i<nd;i++){
-	size_t base_i  = rem % RegionSize[i]; rem /= RegionSize[i];
-	from_coor[i] = base_i + FromLowerLeft[i];
-	to_coor[i] = base_i + ToLowerLeft[i];
-      }
-      
-      int foidx = Fg->oIndex(from_coor);
-      int fiidx = Fg->iIndex(from_coor);
-      int toidx = Tg->oIndex(to_coor);
-      int tiidx = Tg->iIndex(to_coor);
-      int* tt = table + 4*idx;
-      tt[0] = foidx;
-      tt[1] = fiidx;
-      tt[2] = toidx;
-      tt[3] = tiidx;
-    });
-  
-  int* table_d = (int*)acceleratorAllocDevice(tbytes);
-  acceleratorCopyToDevice(table,table_d,tbytes);
-
-  typedef typename vobj::vector_type vector_type;
-  typedef typename vobj::scalar_type scalar_type;
-
-  autoView(from_v,From,AcceleratorRead);
-  autoView(to_v,To,AcceleratorWrite);
-  
-  accelerator_for(idx,nsite,1,{
-      static const int words=sizeof(vobj)/sizeof(vector_type);
-      int* tt = table_d + 4*idx;
-      int from_oidx = *tt++;
-      int from_lane = *tt++;
-      int to_oidx = *tt++;
-      int to_lane = *tt;
-
-      const vector_type* from = (const vector_type *)&from_v[from_oidx];
-      vector_type* to = (vector_type *)&to_v[to_oidx];
-      
-      scalar_type stmp;
-      for(int w=0;w<words;w++){
-	stmp = getlane(from[w], from_lane);
-	putlane(to[w], stmp, to_lane);
-      }
-    });
-  
-  acceleratorFreeDevice(table_d);    
-  free(table);
-  
-
-#else  
   Coordinate ldf = Fg->_ldimensions;
   Coordinate rdf = Fg->_rdimensions;
   Coordinate isf = Fg->_istride;
@@ -765,9 +707,9 @@ void localCopyRegion(const Lattice<vobj> &From,Lattice<vobj> & To,Coordinate Fro
   Coordinate ist = Tg->_istride;
   Coordinate ost = Tg->_ostride;
 
-  autoView( t_v , To, CpuWrite);
+  autoView( t_v , To, AcceleratorWrite);
-  autoView( f_v , From, CpuRead);
+  autoView( f_v , From, AcceleratorRead);
-  thread_for(idx,Fg->lSites(),{
+  accelerator_for(idx,Fg->lSites(),1,{
     sobj s;
     Coordinate Fcoor(nd);
     Coordinate Tcoor(nd);
@@ -780,24 +722,17 @@ void localCopyRegion(const Lattice<vobj> &From,Lattice<vobj> & To,Coordinate Fro
       Tcoor[d] = ToLowerLeft[d]+ Fcoor[d]-FromLowerLeft[d];
     }
     if (in_region) {
-#if 0      
+      Integer idx_f = 0; for(int d=0;d<nd;d++) idx_f+=isf[d]*(Fcoor[d]/rdf[d]);
-      Integer idx_f = 0; for(int d=0;d<nd;d++) idx_f+=isf[d]*(Fcoor[d]/rdf[d]); // inner index from
+      Integer idx_t = 0; for(int d=0;d<nd;d++) idx_t+=ist[d]*(Tcoor[d]/rdt[d]);
-      Integer idx_t = 0; for(int d=0;d<nd;d++) idx_t+=ist[d]*(Tcoor[d]/rdt[d]); // inner index to
+      Integer odx_f = 0; for(int d=0;d<nd;d++) odx_f+=osf[d]*(Fcoor[d]%rdf[d]);
-      Integer odx_f = 0; for(int d=0;d<nd;d++) odx_f+=osf[d]*(Fcoor[d]%rdf[d]); // outer index from
+      Integer odx_t = 0; for(int d=0;d<nd;d++) odx_t+=ost[d]*(Tcoor[d]%rdt[d]);
-      Integer odx_t = 0; for(int d=0;d<nd;d++) odx_t+=ost[d]*(Tcoor[d]%rdt[d]); // outer index to
+      vector_type * fp = (vector_type *)&f_v[odx_f];
-      scalar_type * fp = (scalar_type *)&f_v[odx_f];
+      vector_type * tp = (vector_type *)&t_v[odx_t];
-      scalar_type * tp = (scalar_type *)&t_v[odx_t];
       for(int w=0;w<words;w++){
 	tp[w].putlane(fp[w].getlane(idx_f),idx_t);
       }
-#else
-    peekLocalSite(s,f_v,Fcoor);
-    pokeLocalSite(s,t_v,Tcoor);
-#endif
     }
   });
-
-#endif
 }
 
 
@@ -890,8 +825,6 @@ void ExtractSlice(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int slic
 }
 
 
-//Insert subvolume orthogonal to direction 'orthog' with slice index 'slice_lo' from 'lowDim' onto slice index 'slice_hi' of higherDim
-//The local dimensions of both 'lowDim' and 'higherDim' orthogonal to 'orthog' should be the same
 template<class vobj>
 void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
 {
@@ -908,70 +841,11 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
 
   for(int d=0;d<nh;d++){
     if ( d!=orthog ) {
-      assert(lg->_processors[d]  == hg->_processors[d]);
+    assert(lg->_processors[d]  == hg->_processors[d]);
-      assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
+    assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
-    }
+  }
   }
 
-#if 1
-  size_t nsite = lg->lSites()/lg->LocalDimensions()[orthog];
-  size_t tbytes = 4*nsite*sizeof(int);
-  int *table = (int*)malloc(tbytes);
-  
-  thread_for(idx,nsite,{
-    Coordinate lcoor(nl);
-    Coordinate hcoor(nh);
-    lcoor[orthog] = slice_lo;
-    hcoor[orthog] = slice_hi;
-    size_t rem = idx;
-    for(int mu=0;mu<nl;mu++){
-      if(mu != orthog){
-	int xmu = rem % lg->LocalDimensions()[mu];  rem /= lg->LocalDimensions()[mu];
-	lcoor[mu] = hcoor[mu] = xmu;
-      }
-    }
-    int loidx = lg->oIndex(lcoor);
-    int liidx = lg->iIndex(lcoor);
-    int hoidx = hg->oIndex(hcoor);
-    int hiidx = hg->iIndex(hcoor);
-    int* tt = table + 4*idx;
-    tt[0] = loidx;
-    tt[1] = liidx;
-    tt[2] = hoidx;
-    tt[3] = hiidx;
-    });
-   
-  int* table_d = (int*)acceleratorAllocDevice(tbytes);
-  acceleratorCopyToDevice(table,table_d,tbytes);
-
-  typedef typename vobj::vector_type vector_type;
-  typedef typename vobj::scalar_type scalar_type;
-
-  autoView(lowDim_v,lowDim,AcceleratorRead);
-  autoView(higherDim_v,higherDim,AcceleratorWrite);
-  
-  accelerator_for(idx,nsite,1,{
-      static const int words=sizeof(vobj)/sizeof(vector_type);
-      int* tt = table_d + 4*idx;
-      int from_oidx = *tt++;
-      int from_lane = *tt++;
-      int to_oidx = *tt++;
-      int to_lane = *tt;
-
-      const vector_type* from = (const vector_type *)&lowDim_v[from_oidx];
-      vector_type* to = (vector_type *)&higherDim_v[to_oidx];
-      
-      scalar_type stmp;
-      for(int w=0;w<words;w++){
-	stmp = getlane(from[w], from_lane);
-	putlane(to[w], stmp, to_lane);
-      }
-    });
-  
-  acceleratorFreeDevice(table_d);    
-  free(table);
-  
-#else
   // the above should guarantee that the operations are local
   autoView(lowDimv,lowDim,CpuRead);
   autoView(higherDimv,higherDim,CpuWrite);
@@ -987,7 +861,6 @@ void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int
       pokeLocalSite(s,higherDimv,hcoor);
     }
   });
-#endif
 }
 
 

Grid/lattice/Lattice_view.h

@@ -45,7 +45,6 @@ public:
   };
   // Host only
   GridBase * getGrid(void) const { return _grid; };
-  vobj* getHostPointer(void) const { return _odata; };
 };
 
 /////////////////////////////////////////////////////////////////////////////////////////

Grid/lattice/PaddedCell.h

@@ -1,174 +0,0 @@
-/*************************************************************************************
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/lattice/PaddedCell.h
-
-    Copyright (C) 2019
-
-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
-
-#include<Grid/cshift/Cshift.h>
-
-NAMESPACE_BEGIN(Grid);
-
-//Allow the user to specify how the C-shift is performed, e.g. to respect the appropriate boundary conditions
-template<typename vobj>
-struct CshiftImplBase{
-  virtual Lattice<vobj> Cshift(const Lattice<vobj> &in, int dir, int shift) const = 0;
-  virtual ~CshiftImplBase(){}
-};
-template<typename vobj>
-struct CshiftImplDefault: public CshiftImplBase<vobj>{
-  Lattice<vobj> Cshift(const Lattice<vobj> &in, int dir, int shift) const override{ return Grid::Cshift(in,dir,shift); }
-};
-template<typename Gimpl>
-struct CshiftImplGauge: public CshiftImplBase<typename Gimpl::GaugeLinkField::vector_object>{
-  typename Gimpl::GaugeLinkField Cshift(const typename Gimpl::GaugeLinkField &in, int dir, int shift) const override{ return Gimpl::CshiftLink(in,dir,shift); }
-};  
-
-class PaddedCell {
-public:
-  GridCartesian * unpadded_grid;
-  int dims;
-  int depth;
-  std::vector<GridCartesian *> grids;
-
-  ~PaddedCell()
-  {
-    DeleteGrids();
-  }
-  PaddedCell(int _depth,GridCartesian *_grid)
-  {
-    unpadded_grid = _grid;
-    depth=_depth;
-    dims=_grid->Nd();
-    AllocateGrids();
-    Coordinate local     =unpadded_grid->LocalDimensions();
-    for(int d=0;d<dims;d++){
-      assert(local[d]>=depth);
-    }
-  }
-  void DeleteGrids(void)
-  {
-    for(int d=0;d<grids.size();d++){
-      delete grids[d];
-    }
-    grids.resize(0);
-  };
-  void AllocateGrids(void)
-  {
-    Coordinate local     =unpadded_grid->LocalDimensions();
-    Coordinate simd      =unpadded_grid->_simd_layout;
-    Coordinate processors=unpadded_grid->_processors;
-    Coordinate plocal    =unpadded_grid->LocalDimensions();
-    Coordinate global(dims);
-
-    // expand up one dim at a time
-    for(int d=0;d<dims;d++){
-
-      plocal[d] += 2*depth; 
-
-      for(int d=0;d<dims;d++){
-	global[d] = plocal[d]*processors[d];
-      }
-
-      grids.push_back(new GridCartesian(global,simd,processors));
-    }
-  };
-  template<class vobj>
-  inline Lattice<vobj> Extract(const Lattice<vobj> &in) const
-  {
-    Lattice<vobj> out(unpadded_grid);
-
-    Coordinate local     =unpadded_grid->LocalDimensions();
-    Coordinate fll(dims,depth); // depends on the MPI spread
-    Coordinate tll(dims,0); // depends on the MPI spread
-    localCopyRegion(in,out,fll,tll,local);
-    return out;
-  }
-  template<class vobj>
-  inline Lattice<vobj> Exchange(const Lattice<vobj> &in, const CshiftImplBase<vobj> &cshift = CshiftImplDefault<vobj>()) const
-  {
-    GridBase *old_grid = in.Grid();
-    int dims = old_grid->Nd();
-    Lattice<vobj> tmp = in;
-    for(int d=0;d<dims;d++){
-      tmp = Expand(d,tmp,cshift); // rvalue && assignment
-    }
-    return tmp;
-  }
-  // expand up one dim at a time
-  template<class vobj>
-  inline Lattice<vobj> Expand(int dim, const Lattice<vobj> &in, const CshiftImplBase<vobj> &cshift = CshiftImplDefault<vobj>()) const
-  {
-    GridBase *old_grid = in.Grid();
-    GridCartesian *new_grid = grids[dim];//These are new grids
-    Lattice<vobj>  padded(new_grid);
-    Lattice<vobj> shifted(old_grid);    
-    Coordinate local     =old_grid->LocalDimensions();
-    Coordinate plocal    =new_grid->LocalDimensions();
-    if(dim==0) conformable(old_grid,unpadded_grid);
-    else       conformable(old_grid,grids[dim-1]);
-
-    std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
-
-    double tins=0, tshift=0;
-    
-    // Middle bit
-    double t = usecond();
-    for(int x=0;x<local[dim];x++){
-      InsertSliceLocal(in,padded,x,depth+x,dim);
-    }
-    tins += usecond() - t;
-    
-    // High bit
-    t = usecond();
-    shifted = cshift.Cshift(in,dim,depth);
-    tshift += usecond() - t;
-
-    t=usecond();
-    for(int x=0;x<depth;x++){
-      InsertSliceLocal(shifted,padded,local[dim]-depth+x,depth+local[dim]+x,dim);
-    }
-    tins += usecond() - t;
-    
-    // Low bit
-    t = usecond();
-    shifted = cshift.Cshift(in,dim,-depth);
-    tshift += usecond() - t;
-    
-    t = usecond();
-    for(int x=0;x<depth;x++){
-      InsertSliceLocal(shifted,padded,x,x,dim);
-    }
-    tins += usecond() - t;
-
-    std::cout << GridLogPerformance << "PaddedCell::Expand timings: cshift:" << tshift/1000 << "ms, insert-slice:" << tins/1000 << "ms" << std::endl;
-    
-    return padded;
-  }
-
-};
- 
-
-NAMESPACE_END(Grid);
-

Grid/log/Log.h

@@ -179,11 +179,11 @@ extern GridLogger GridLogSolver;
 extern GridLogger GridLogError;
 extern GridLogger GridLogWarning;
 extern GridLogger GridLogMessage;
-extern GridLogger GridLogDebug;
+extern GridLogger GridLogDebug  ;
 extern GridLogger GridLogPerformance;
 extern GridLogger GridLogDslash;
-extern GridLogger GridLogIterative;
+extern GridLogger GridLogIterative  ;
-extern GridLogger GridLogIntegrator;
+extern GridLogger GridLogIntegrator  ;
 extern GridLogger GridLogHMC;
 extern GridLogger GridLogMemory;
 extern GridLogger GridLogTracing;
@@ -191,41 +191,6 @@ extern Colours    GridLogColours;
 
 std::string demangle(const char* name) ;
 
-template<typename... Args>
-inline std::string sjoin(Args&&... args) noexcept {
-    std::ostringstream msg;
-    (msg << ... << args);
-    return msg.str();
-}
-
-/*!  @brief make log messages work like python print */
-template <typename... Args>
-inline void Grid_log(Args&&... args) {
-    std::string msg = sjoin(std::forward<Args>(args)...);
-    std::cout << GridLogMessage << msg << std::endl;
-}
-
-/*!  @brief make warning messages work like python print */
-template <typename... Args>
-inline void Grid_warn(Args&&... args) {
-    std::string msg = sjoin(std::forward<Args>(args)...);
-    std::cout << "\033[33m" << GridLogWarning << msg << "\033[0m" << std::endl;
-}
-
-/*!  @brief make error messages work like python print */
-template <typename... Args>
-inline void Grid_error(Args&&... args) {
-    std::string msg = sjoin(std::forward<Args>(args)...);
-    std::cout << "\033[31m" << GridLogError << msg << "\033[0m" << std::endl;
-}
-
-/*!  @brief make pass messages work like python print */
-template <typename... Args>
-inline void Grid_pass(Args&&... args) {
-    std::string msg = sjoin(std::forward<Args>(args)...);
-    std::cout << "\033[32m" << GridLogMessage << msg << "\033[0m" << std::endl;
-}
-
 #define _NBACKTRACE (256)
 extern void * Grid_backtrace_buffer[_NBACKTRACE];
 

Grid/perfmon/Tracing.h

@@ -34,7 +34,7 @@ class GridTracer {
 };
 inline void tracePush(const char *name) { roctxRangePushA(name); }
 inline void tracePop(const char *name) { roctxRangePop(); }
-inline int  traceStart(const char *name) { return roctxRangeStart(name); }
+inline int  traceStart(const char *name) { roctxRangeStart(name); }
 inline void traceStop(int ID) { roctxRangeStop(ID); }
 #endif
 

Grid/qcd/QCD.h

@@ -104,7 +104,6 @@ template<typename vtype> using iSpinMatrix                = iScalar<iMatrix<iSca
 template<typename vtype> using iColourMatrix              = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
 template<typename vtype> using iSpinColourMatrix          = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
 template<typename vtype> using iLorentzColourMatrix       = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
-template<typename vtype> using iLorentzComplex            = iVector<iScalar<iScalar<vtype> >, Nd > ;
 template<typename vtype> using iDoubleStoredColourMatrix  = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
 template<typename vtype> using iSpinVector                = iScalar<iVector<iScalar<vtype>, Ns> >;
 template<typename vtype> using iColourVector              = iScalar<iScalar<iVector<vtype, Nc> > >;
@@ -179,15 +178,6 @@ typedef iLorentzColourMatrix<vComplexF>  vLorentzColourMatrixF;
 typedef iLorentzColourMatrix<vComplexD>  vLorentzColourMatrixD;
 typedef iLorentzColourMatrix<vComplexD2> vLorentzColourMatrixD2;
 
-// LorentzComplex
-typedef iLorentzComplex<Complex  > LorentzComplex;
-typedef iLorentzComplex<ComplexF > LorentzComplexF;
-typedef iLorentzComplex<ComplexD > LorentzComplexD;
-
-typedef iLorentzComplex<vComplex > vLorentzComplex;
-typedef iLorentzComplex<vComplexF> vLorentzComplexF;
-typedef iLorentzComplex<vComplexD> vLorentzComplexD;
-
 // DoubleStored gauge field
 typedef iDoubleStoredColourMatrix<Complex  > DoubleStoredColourMatrix;
 typedef iDoubleStoredColourMatrix<ComplexF > DoubleStoredColourMatrixF;
@@ -317,10 +307,6 @@ typedef Lattice<vLorentzColourMatrixF>  LatticeLorentzColourMatrixF;
 typedef Lattice<vLorentzColourMatrixD>  LatticeLorentzColourMatrixD;
 typedef Lattice<vLorentzColourMatrixD2> LatticeLorentzColourMatrixD2;
 
-typedef Lattice<vLorentzComplex>  LatticeLorentzComplex;
-typedef Lattice<vLorentzComplexF> LatticeLorentzComplexF;
-typedef Lattice<vLorentzComplexD> LatticeLorentzComplexD;
-
 // DoubleStored gauge field
 typedef Lattice<vDoubleStoredColourMatrix>   LatticeDoubleStoredColourMatrix;
 typedef Lattice<vDoubleStoredColourMatrixF>  LatticeDoubleStoredColourMatrixF;

Grid/qcd/action/ActionBase.h

@@ -34,24 +34,10 @@ directory
 
 NAMESPACE_BEGIN(Grid);
 
-///////////////////////////////////
-// Smart configuration base class
-///////////////////////////////////
-template< class Field >
-class ConfigurationBase
-{
-public:
-  ConfigurationBase() {}
-  virtual ~ConfigurationBase() {}
-  virtual void set_Field(Field& U) =0;
-  virtual void smeared_force(Field&) = 0;
-  virtual Field& get_SmearedU() =0;
-  virtual Field &get_U(bool smeared = false) = 0;
-};
-
 template <class GaugeField >
 class Action 
 {
+
 public:
   bool is_smeared = false;
   RealD deriv_norm_sum;
@@ -91,60 +77,16 @@ public:
   void refresh_timer_stop(void)  { refresh_us+=usecond(); }
   void S_timer_start(void)       { S_us-=usecond(); }
   void S_timer_stop(void)        { S_us+=usecond(); }
-  /////////////////////////////
   // Heatbath?
-  /////////////////////////////
   virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) = 0; // refresh pseudofermions
   virtual RealD S(const GaugeField& U) = 0;                             // evaluate the action
   virtual RealD Sinitial(const GaugeField& U) { return this->S(U); } ;  // if the refresh computes the action, can cache it. Alternately refreshAndAction() ?
   virtual void deriv(const GaugeField& U, GaugeField& dSdU) = 0;        // evaluate the action derivative
-
-  /////////////////////////////////////////////////////////////
-  // virtual smeared interface through configuration container
-  /////////////////////////////////////////////////////////////
-  virtual void refresh(ConfigurationBase<GaugeField> & U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
-  {
-    refresh(U.get_U(is_smeared),sRNG,pRNG);
-  }
-  virtual RealD S(ConfigurationBase<GaugeField>& U)
-  {
-    return S(U.get_U(is_smeared));
-  }
-  virtual RealD Sinitial(ConfigurationBase<GaugeField>& U) 
-  {
-    return Sinitial(U.get_U(is_smeared));
-  }
-  virtual void deriv(ConfigurationBase<GaugeField>& U, GaugeField& dSdU)
-  {
-    deriv(U.get_U(is_smeared),dSdU); 
-    if ( is_smeared ) {
-      U.smeared_force(dSdU);
-    }
-  }
-  ///////////////////////////////
-  // Logging
-  ///////////////////////////////
   virtual std::string action_name()    = 0;                             // return the action name
   virtual std::string LogParameters()  = 0;                             // prints action parameters
   virtual ~Action(){}
 };
 
-template <class GaugeField >
-class EmptyAction : public Action <GaugeField>
-{
-  virtual void refresh(const GaugeField& U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) { assert(0);}; // refresh pseudofermions
-  virtual RealD S(const GaugeField& U) { return 0.0;};                             // evaluate the action
-  virtual void deriv(const GaugeField& U, GaugeField& dSdU) { assert(0); };        // evaluate the action derivative
-
-  ///////////////////////////////
-  // Logging
-  ///////////////////////////////
-  virtual std::string action_name()    { return std::string("Level Force Log"); };
-  virtual std::string LogParameters()  { return std::string("No parameters");};
-};
-
-
-
 NAMESPACE_END(Grid);
 
 #endif // ACTION_BASE_H

Grid/qcd/action/ActionCore.h

@@ -30,8 +30,6 @@ directory
 #ifndef QCD_ACTION_CORE
 #define QCD_ACTION_CORE
 
-#include <Grid/qcd/action/gauge/GaugeImplementations.h>
-
 #include <Grid/qcd/action/ActionBase.h>
 NAMESPACE_CHECK(ActionBase);
 #include <Grid/qcd/action/ActionSet.h>

Grid/qcd/action/fermion/Fermion.h

@@ -126,16 +126,6 @@ typedef WilsonFermion<WilsonTwoIndexSymmetricImplD> WilsonTwoIndexSymmetricFermi
 typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplF> WilsonTwoIndexAntiSymmetricFermionF;
 typedef WilsonFermion<WilsonTwoIndexAntiSymmetricImplD> WilsonTwoIndexAntiSymmetricFermionD;
 
-// Sp(2n)
-typedef WilsonFermion<SpWilsonImplF> SpWilsonFermionF;
-typedef WilsonFermion<SpWilsonImplD> SpWilsonFermionD;
-
-typedef WilsonFermion<SpWilsonTwoIndexAntiSymmetricImplF> SpWilsonTwoIndexAntiSymmetricFermionF;
-typedef WilsonFermion<SpWilsonTwoIndexAntiSymmetricImplD> SpWilsonTwoIndexAntiSymmetricFermionD;
-
-typedef WilsonFermion<SpWilsonTwoIndexSymmetricImplF> SpWilsonTwoIndexSymmetricFermionF;
-typedef WilsonFermion<SpWilsonTwoIndexSymmetricImplD> SpWilsonTwoIndexSymmetricFermionD;
-
 // Twisted mass fermion
 typedef WilsonTMFermion<WilsonImplD2> WilsonTMFermionD2;
 typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;

Grid/qcd/action/fermion/WilsonImpl.h

@@ -261,22 +261,6 @@ typedef WilsonImpl<vComplex,  TwoIndexAntiSymmetricRepresentation, CoeffReal > W
 typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF;  // Float
 typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD;  // Double
 
-//sp 2n
-
-typedef WilsonImpl<vComplex,  SpFundamentalRepresentation, CoeffReal > SpWilsonImplR;  // Real.. whichever prec
-typedef WilsonImpl<vComplexF, SpFundamentalRepresentation, CoeffReal > SpWilsonImplF;  // Float
-typedef WilsonImpl<vComplexD, SpFundamentalRepresentation, CoeffReal > SpWilsonImplD;  // Double
-
-typedef WilsonImpl<vComplex,  SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplR;  // Real.. whichever prec
-typedef WilsonImpl<vComplexF, SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplF;  // Float
-typedef WilsonImpl<vComplexD, SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplD;  // Double
-
-typedef WilsonImpl<vComplex,  SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplR;  // Real.. whichever prec
-typedef WilsonImpl<vComplexF, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplF;  // Float
-typedef WilsonImpl<vComplexD, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplD;  // Double
-
-typedef WilsonImpl<vComplex,  SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonAdjImplR;  // Real.. whichever prec    // adj = 2indx symmetric for Sp(2N)
-typedef WilsonImpl<vComplexF, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonAdjImplF;  // Float     // adj = 2indx symmetric for Sp(2N)
-typedef WilsonImpl<vComplexD, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonAdjImplD;  // Double    // adj = 2indx symmetric for Sp(2N)
 
 NAMESPACE_END(Grid);
+

Grid/qcd/action/fermion/WilsonTMFermion.h

@@ -63,9 +63,7 @@ public:
   virtual void MooeeDag(const FermionField &in, FermionField &out) ;
   virtual void MooeeInv(const FermionField &in, FermionField &out) ;
   virtual void MooeeInvDag(const FermionField &in, FermionField &out) ;
-  virtual void M(const FermionField &in, FermionField &out) ;
+
-  virtual void Mdag(const FermionField &in, FermionField &out) ;
-  
 private:
   RealD mu; // TwistedMass parameter
 

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/WilsonFermion5DImplementation.h

@@ -196,7 +196,6 @@ void WilsonFermion5D<Impl>::DhopDir(const FermionField &in, FermionField &out,in
   
   uint64_t Nsite = Umu.Grid()->oSites();
   Kernels::DhopDirKernel(Stencil,Umu,Stencil.CommBuf(),Ls,Nsite,in,out,dirdisp,gamma);
-
 };
 template<class Impl>
 void WilsonFermion5D<Impl>::DhopDirAll(const FermionField &in, std::vector<FermionField> &out)
@@ -247,10 +246,14 @@ void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
 
     Kernels::DhopDirKernel(st, U, st.CommBuf(), Ls, Usites, B, Btilde, mu,gamma);
 
+    std::cout << " InsertForce Btilde "<< norm2(Btilde)<<std::endl;
+
     ////////////////////////////
     // spin trace outer product
     ////////////////////////////
     Impl::InsertForce5D(mat, Btilde, Atilde, mu);
+
+    std::cout << " InsertForce "<< norm2(mat)<<std::endl;
   }
 }
 

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

@@ -423,6 +423,7 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
 #define KERNEL_CALL(A) KERNEL_CALLNB(A); accelerator_barrier();
 
 #define KERNEL_CALL_EXT(A)						\
+  const uint64_t    NN = Nsite*Ls;					\
   const uint64_t    sz = st.surface_list.size();			\
   auto ptr = &st.surface_list[0];					\
   accelerator_forNB( ss, sz, Simd::Nsimd(), {				\
@@ -462,7 +463,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 +496,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/action/fermion/implementation/WilsonTMFermionImplementation.h

@@ -93,25 +93,5 @@ void WilsonTMFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &ou
   RealD b    = tm /sq;
   axpibg5x(out,in,a,b);
 }
-template<class Impl>
-void WilsonTMFermion<Impl>::M(const FermionField &in, FermionField &out) {
-  out.Checkerboard() = in.Checkerboard();
-  this->Dhop(in, out, DaggerNo);
-  FermionField tmp(out.Grid());
-  RealD a = 4.0+this->mass;
-  RealD b = this->mu;
-  axpibg5x(tmp,in,a,b);
-  axpy(out, 1.0, tmp, out);
-}
-template<class Impl>
-void WilsonTMFermion<Impl>::Mdag(const FermionField &in, FermionField &out) {
-  out.Checkerboard() = in.Checkerboard();
-  this->Dhop(in, out, DaggerYes);
-  FermionField tmp(out.Grid());
-  RealD a = 4.0+this->mass;
-  RealD b = -this->mu;
-  axpibg5x(tmp,in,a,b);
-  axpy(out, 1.0, tmp, out);
-}
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonCloverFermionInstantiationSpWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonFermionInstantiationSpWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonKernelsInstantiationSpWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplD/WilsonTMFermionInstantiationSpWilsonImplD.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplD/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonImplD

Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonCloverFermionInstantiationSpWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonFermionInstantiationSpWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonKernelsInstantiationSpWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplF/WilsonTMFermionInstantiationSpWilsonImplF.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonImplF/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonImplF

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonCloverFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonTMFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonTwoIndexAntiSymmetricImplD

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonCloverFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonTMFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonTwoIndexAntiSymmetricImplF

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonCloverFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationSpWilsonTwoIndexSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonTMFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonTwoIndexSymmetricImplD

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonCloverFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonCloverFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationSpWilsonTwoIndexSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonKernelsInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonTMFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc

@@ -1 +0,0 @@
-../WilsonTMFermionInstantiation.cc.master

Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/impl.h

@@ -1 +0,0 @@
-#define IMPLEMENTATION SpWilsonTwoIndexSymmetricImplF

Grid/qcd/action/fermion/instantiation/generate_instantiations.sh

@@ -10,18 +10,12 @@ WILSON_IMPL_LIST=" \
 	   WilsonImplF \
 	   WilsonImplD \
 	   WilsonImplD2 \
-	   SpWilsonImplF \
-	   SpWilsonImplD \
 	   WilsonAdjImplF \
 	   WilsonAdjImplD \
 	   WilsonTwoIndexSymmetricImplF \
 	   WilsonTwoIndexSymmetricImplD \
 	   WilsonTwoIndexAntiSymmetricImplF \
 	   WilsonTwoIndexAntiSymmetricImplD \
-	   SpWilsonTwoIndexAntiSymmetricImplF \
-	   SpWilsonTwoIndexAntiSymmetricImplD \
-	   SpWilsonTwoIndexSymmetricImplF \
-	   SpWilsonTwoIndexSymmetricImplD \
 	   GparityWilsonImplF \
 	   GparityWilsonImplD "
 

Grid/qcd/action/gauge/Gauge.h

@@ -39,9 +39,6 @@ NAMESPACE_BEGIN(Grid);
 typedef WilsonGaugeAction<PeriodicGimplR>          WilsonGaugeActionR;
 typedef WilsonGaugeAction<PeriodicGimplF>          WilsonGaugeActionF;
 typedef WilsonGaugeAction<PeriodicGimplD>          WilsonGaugeActionD;
-typedef WilsonGaugeAction<SpPeriodicGimplR>        SpWilsonGaugeActionR;
-typedef WilsonGaugeAction<SpPeriodicGimplF>        SpWilsonGaugeActionF;
-typedef WilsonGaugeAction<SpPeriodicGimplD>        SpWilsonGaugeActionD;
 typedef PlaqPlusRectangleAction<PeriodicGimplR>    PlaqPlusRectangleActionR;
 typedef PlaqPlusRectangleAction<PeriodicGimplF>    PlaqPlusRectangleActionF;
 typedef PlaqPlusRectangleAction<PeriodicGimplD>    PlaqPlusRectangleActionD;

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -61,7 +61,7 @@ NAMESPACE_BEGIN(Grid);
   typedef typename Impl::Field Field;
 
 // hardcodes the exponential approximation in the template
-template <class S, int Nrepresentation = Nc, int Nexp = 12, class Group = SU<Nc> > class GaugeImplTypes {
+template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
 public:
   typedef S Simd;
   typedef typename Simd::scalar_type scalar_type;
@@ -78,6 +78,8 @@ public:
   typedef Lattice<SiteLink>    LinkField; 
   typedef Lattice<SiteField>   Field;
 
+  typedef SU<Nrepresentation> Group;
+
   // Guido: we can probably separate the types from the HMC functions
   // this will create 2 kind of implementations
   // probably confusing the users
@@ -117,7 +119,6 @@ public:
     //
     LinkField Pmu(P.Grid());
     Pmu = Zero();
-
     for (int mu = 0; mu < Nd; mu++) {
       Group::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
       RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
@@ -125,12 +126,8 @@ public:
       PokeIndex<LorentzIndex>(P, Pmu, mu);
     }
   }
-    
+
-  static inline Field projectForce(Field &P) {
+  static inline Field projectForce(Field &P) { return Ta(P); }
-      Field ret(P.Grid());
-      Group::taProj(P, ret);
-      return ret;
-    }
 
   static inline void update_field(Field& P, Field& U, double ep){
     //static std::chrono::duration<double> diff;
@@ -140,15 +137,14 @@ public:
     autoView(P_v,P,AcceleratorRead);
     accelerator_for(ss, P.Grid()->oSites(),1,{
       for (int mu = 0; mu < Nd; mu++) {
-          U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu);
+        U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
-          U_v[ss](mu) = Group::ProjectOnGeneralGroup(U_v[ss](mu));
       }
     });
    //auto end = std::chrono::high_resolution_clock::now();
    // diff += end - start;
    // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n";
   }
-    
+
   static inline RealD FieldSquareNorm(Field& U){
     LatticeComplex Hloc(U.Grid());
     Hloc = Zero();
@@ -161,7 +157,7 @@ public:
   }
 
   static inline void Project(Field &U) {
-    Group::ProjectOnSpecialGroup(U);
+    ProjectSUn(U);
   }
 
   static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
@@ -175,7 +171,6 @@ public:
   static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
     Group::ColdConfiguration(pRNG, U);
   }
-
 };
 
 
@@ -183,17 +178,10 @@ typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
 typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
 typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
 
-typedef GaugeImplTypes<vComplex, Nc, 12, Sp<Nc> > SpGimplTypesR;
-typedef GaugeImplTypes<vComplexF, Nc, 12, Sp<Nc> > SpGimplTypesF;
-typedef GaugeImplTypes<vComplexD, Nc, 12, Sp<Nc> > SpGimplTypesD;
-
 typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
 typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
 typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
 
-
-
-
 NAMESPACE_END(Grid);
 
 #endif // GRID_GAUGE_IMPL_TYPES_H

Grid/qcd/action/gauge/GaugeImplementations.h

@@ -176,7 +176,7 @@ public:
       return PeriodicBC::CshiftLink(Link,mu,shift);
   }
 
-  static inline void       setDirections(const std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
+  static inline void       setDirections(std::vector<int> &conjDirs) { _conjDirs=conjDirs; }
   static inline std::vector<int> getDirections(void) { return _conjDirs; }
   static inline bool isPeriodicGaugeField(void) { return false; }
 };
@@ -193,11 +193,6 @@ typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever pre
 typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
 typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
 
-typedef PeriodicGaugeImpl<SpGimplTypesR> SpPeriodicGimplR; // Real.. whichever prec
-typedef PeriodicGaugeImpl<SpGimplTypesF> SpPeriodicGimplF; // Float
-typedef PeriodicGaugeImpl<SpGimplTypesD> SpPeriodicGimplD; // Double
-
-
 NAMESPACE_END(Grid);
 
 #endif

Grid/qcd/action/gauge/PlaqPlusRectangleAction.h

@@ -43,7 +43,7 @@ public:
 private:
   RealD c_plaq;
   RealD c_rect;
-  typename WilsonLoops<Gimpl>::StapleAndRectStapleAllWorkspace workspace;
+
 public:
   PlaqPlusRectangleAction(RealD b,RealD c): c_plaq(b),c_rect(c){};
 
@@ -79,18 +79,27 @@ public:
     GridBase *grid = Umu.Grid();
 
     std::vector<GaugeLinkField> U (Nd,grid);
+    std::vector<GaugeLinkField> U2(Nd,grid);
+
     for(int mu=0;mu<Nd;mu++){
       U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+      WilsonLoops<Gimpl>::RectStapleDouble(U2[mu],U[mu],mu);
     }
-    std::vector<GaugeLinkField> RectStaple(Nd,grid), Staple(Nd,grid);
-    WilsonLoops<Gimpl>::StapleAndRectStapleAll(Staple, RectStaple, U, workspace);
 
     GaugeLinkField dSdU_mu(grid);
     GaugeLinkField staple(grid);
 
     for (int mu=0; mu < Nd; mu++){
-      dSdU_mu = Ta(U[mu]*Staple[mu])*factor_p;
+
-      dSdU_mu = dSdU_mu + Ta(U[mu]*RectStaple[mu])*factor_r;
+      // Staple in direction mu
+
+      WilsonLoops<Gimpl>::Staple(staple,Umu,mu);
+
+      dSdU_mu = Ta(U[mu]*staple)*factor_p;
+
+      WilsonLoops<Gimpl>::RectStaple(Umu,staple,U2,U,mu);
+
+      dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r;
 	  
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }

Grid/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h

@@ -119,13 +119,19 @@ public:
     //  X^dag Der_oe MeeInv Meo Y
     // Use Mooee as nontrivial but gauge field indept
     this->_Mat.MeooeDag   (V,tmp1);      // odd->even -- implicit -0.5 factor to be applied
+    std::cout << " tmp 1" << norm2(tmp1)<<std::endl;
     this->_Mat.MooeeInvDag(tmp1,tmp2);   // even->even 
+    std::cout << " tmp 1" << norm2(tmp2)<<std::endl;
     this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes);
+    std::cout << " ForceO " << norm2(ForceO)<<std::endl;
           
     //  Accumulate X^dag M_oe MeeInv Der_eo Y
     this->_Mat.Meooe   (U,tmp1);    // even->odd -- implicit -0.5 factor to be applied
+    std::cout << " tmp 1" << norm2(tmp1)<<std::endl;
     this->_Mat.MooeeInv(tmp1,tmp2); // even->even 
+    std::cout << " tmp 2" << norm2(tmp2)<<std::endl;
     this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes);
+    std::cout << " ForceE " << norm2(ForceE)<<std::endl;
 
     assert(ForceE.Checkerboard()==Even);
     assert(ForceO.Checkerboard()==Odd);

Grid/qcd/hmc/GenericHMCrunner.h

@@ -225,18 +225,6 @@ template <class RepresentationsPolicy,
 using GenericHMCRunnerHirep =
 				     HMCWrapperTemplate<PeriodicGimplR, Integrator, RepresentationsPolicy>;
 
-// sp2n
-
-template <template <typename, typename, typename> class Integrator>
-using GenericSpHMCRunner = HMCWrapperTemplate<SpPeriodicGimplR, Integrator>;
-
-template <class RepresentationsPolicy,
-          template <typename, typename, typename> class Integrator>
-using GenericSpHMCRunnerHirep =
-                     HMCWrapperTemplate<SpPeriodicGimplR, Integrator, RepresentationsPolicy>;
-
-
-
 template <class Implementation, class RepresentationsPolicy, 
           template <typename, typename, typename> class Integrator>
 using GenericHMCRunnerTemplate = HMCWrapperTemplate<Implementation, Integrator, RepresentationsPolicy>;

Grid/qcd/hmc/HMC.h

@@ -283,13 +283,12 @@ public:
       std::cout << GridLogHMC << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;
 
       TheIntegrator.print_timer();
-      
+
-      TheIntegrator.Smearer.set_Field(Ucur);
       for (int obs = 0; obs < Observables.size(); obs++) {
       	std::cout << GridLogDebug << "Observables # " << obs << std::endl;
       	std::cout << GridLogDebug << "Observables total " << Observables.size() << std::endl;
       	std::cout << GridLogDebug << "Observables pointer " << Observables[obs] << std::endl;
-        Observables[obs]->TrajectoryComplete(traj + 1, TheIntegrator.Smearer, sRNG, pRNG);
+        Observables[obs]->TrajectoryComplete(traj + 1, Ucur, sRNG, pRNG);
       }
       std::cout << GridLogHMC << ":::::::::::::::::::::::::::::::::::::::::::" << std::endl;
     }

Grid/qcd/hmc/checkpointers/BaseCheckpointer.h

@@ -35,16 +35,13 @@ class CheckpointerParameters : Serializable {
 public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(CheckpointerParameters, 
 				  std::string, config_prefix, 
-				  std::string, smeared_prefix, 
 				  std::string, rng_prefix, 
 				  int, saveInterval, 
-				  bool, saveSmeared, 
 				  std::string, format, );
 
-  CheckpointerParameters(std::string cf = "cfg", std::string sf="cfg_smr" , std::string rn = "rng",
+  CheckpointerParameters(std::string cf = "cfg", std::string rn = "rng",
 			 int savemodulo = 1, const std::string &f = "IEEE64BIG")
     : config_prefix(cf),
-      smeared_prefix(sf),
       rng_prefix(rn),
       saveInterval(savemodulo),
       format(f){};
@@ -64,21 +61,13 @@ template <class Impl>
 class BaseHmcCheckpointer : public HmcObservable<typename Impl::Field> {
 public:
   void build_filenames(int traj, CheckpointerParameters &Params,
-                       std::string &conf_file,
+                       std::string &conf_file, std::string &rng_file) {
-                       std::string &smear_file,
-		       std::string &rng_file) {
     {
       std::ostringstream os;
       os << Params.rng_prefix << "." << traj;
       rng_file = os.str();
     }
 
-    {
-      std::ostringstream os;
-      os << Params.smeared_prefix << "." << traj;
-      smear_file = os.str();
-    }
-
     {
       std::ostringstream os;
       os << Params.config_prefix << "." << traj;
@@ -95,11 +84,6 @@ public:
   }
   virtual void initialize(const CheckpointerParameters &Params) = 0;
 
-  virtual void TrajectoryComplete(int traj,
-                                  typename Impl::Field &U,
-                                  GridSerialRNG &sRNG,
-                                  GridParallelRNG &pRNG) { assert(0); } ; // HMC should pass the smart config with smeared and unsmeared
-  
   virtual void CheckpointRestore(int traj, typename Impl::Field &U,
                                  GridSerialRNG &sRNG,
                                  GridParallelRNG &pRNG) = 0;

Grid/qcd/hmc/checkpointers/BinaryCheckpointer.h

@@ -61,14 +61,11 @@ public:
     fout.close();
   }
 
-  void TrajectoryComplete(int traj,
+  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG) {
-			  ConfigurationBase<Field> &SmartConfig,
-			  GridSerialRNG &sRNG, GridParallelRNG &pRNG)
-  {
 
     if ((traj % Params.saveInterval) == 0) {
-      std::string config, rng, smr;
+      std::string config, rng;
-      this->build_filenames(traj, Params, config, smr, rng);
+      this->build_filenames(traj, Params, config, rng);
 
       uint32_t nersc_csum;
       uint32_t scidac_csuma;
@@ -77,15 +74,9 @@ public:
       BinarySimpleUnmunger<sobj_double, sobj> munge;
       truncate(rng);
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
-      std::cout << GridLogMessage << "Written Binary RNG " << rng
-                << " checksum " << std::hex 
-		<< nersc_csum   <<"/"
-		<< scidac_csuma   <<"/"
-		<< scidac_csumb 
-		<< std::dec << std::endl;
-
       truncate(config);
-      BinaryIO::writeLatticeObject<vobj, sobj_double>(SmartConfig.get_U(false), config, munge, 0, Params.format,
+
+      BinaryIO::writeLatticeObject<vobj, sobj_double>(U, config, munge, 0, Params.format,
 						      nersc_csum,scidac_csuma,scidac_csumb);
 
       std::cout << GridLogMessage << "Written Binary Configuration " << config
@@ -94,18 +85,6 @@ public:
 		<< scidac_csuma   <<"/"
 		<< scidac_csumb 
 		<< std::dec << std::endl;
-
-      if ( Params.saveSmeared ) {
-	truncate(smr);
-	BinaryIO::writeLatticeObject<vobj, sobj_double>(SmartConfig.get_U(true), smr, munge, 0, Params.format,
-							nersc_csum,scidac_csuma,scidac_csumb);
-	std::cout << GridLogMessage << "Written Binary Smeared Configuration " << smr
-                << " checksum " << std::hex 
-		<< nersc_csum   <<"/"
-		<< scidac_csuma   <<"/"
-		<< scidac_csumb 
-		<< std::dec << std::endl;
-      }
     }
 
   };

Grid/qcd/hmc/checkpointers/ILDGCheckpointer.h

@@ -69,27 +69,17 @@ public:
     }
   }
 
-  void TrajectoryComplete(int traj,
+  void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
-			  ConfigurationBase<GaugeField> &SmartConfig,
-			  GridSerialRNG &sRNG,
                           GridParallelRNG &pRNG) {
     if ((traj % Params.saveInterval) == 0) {
-      std::string config, rng, smr;
+      std::string config, rng;
       this->build_filenames(traj, Params, config, rng);
-      GridBase *grid = SmartConfig.get_U(false).Grid();
+      GridBase *grid = U.Grid();
       uint32_t nersc_csum,scidac_csuma,scidac_csumb;
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
-      std::cout << GridLogMessage << "Written BINARY RNG " << rng
-                << " checksum " << std::hex 
-		<< nersc_csum<<"/"
-		<< scidac_csuma<<"/"
-		<< scidac_csumb
-		<< std::dec << std::endl;
-
-      
       IldgWriter _IldgWriter(grid->IsBoss());
       _IldgWriter.open(config);
-      _IldgWriter.writeConfiguration<GaugeStats>(SmartConfig.get_U(false), traj, config, config);
+      _IldgWriter.writeConfiguration<GaugeStats>(U, traj, config, config);
       _IldgWriter.close();
 
       std::cout << GridLogMessage << "Written ILDG Configuration on " << config
@@ -98,21 +88,6 @@ public:
 		<< scidac_csuma<<"/"
 		<< scidac_csumb
 		<< std::dec << std::endl;
-
-      if ( Params.saveSmeared ) { 
-	IldgWriter _IldgWriter(grid->IsBoss());
-	_IldgWriter.open(smr);
-	_IldgWriter.writeConfiguration<GaugeStats>(SmartConfig.get_U(true), traj, config, config);
-	_IldgWriter.close();
-
-	std::cout << GridLogMessage << "Written ILDG Configuration on " << smr
-                << " checksum " << std::hex 
-		<< nersc_csum<<"/"
-		<< scidac_csuma<<"/"
-		<< scidac_csumb
-		<< std::dec << std::endl;
-      }
-
     }
   };
 

Grid/qcd/hmc/checkpointers/NerscCheckpointer.h

@@ -52,29 +52,23 @@ public:
     Params.format = "IEEE64BIG";  // fixed, overwrite any other choice
   }
 
-  virtual void TrajectoryComplete(int traj,
+  void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG,
-                                  ConfigurationBase<GaugeField> &SmartConfig,
+                          GridParallelRNG &pRNG) {
-                                  GridSerialRNG &sRNG,
-                                  GridParallelRNG &pRNG)
-  {
     if ((traj % Params.saveInterval) == 0) {
-      std::string config, rng, smr;
+      std::string config, rng;
-      this->build_filenames(traj, Params, config, smr, rng);
+      this->build_filenames(traj, Params, config, rng);
-      
+
       int precision32 = 1;
       int tworow = 0;
       NerscIO::writeRNGState(sRNG, pRNG, rng);
-      NerscIO::writeConfiguration<GaugeStats>(SmartConfig.get_U(false), config, tworow, precision32);
+      NerscIO::writeConfiguration<GaugeStats>(U, config, tworow, precision32);
-      if ( Params.saveSmeared ) {
-	NerscIO::writeConfiguration<GaugeStats>(SmartConfig.get_U(true), smr, tworow, precision32);
-      }
     }
   };
 
   void CheckpointRestore(int traj, GaugeField &U, GridSerialRNG &sRNG,
                          GridParallelRNG &pRNG) {
-    std::string config, rng, smr;
+    std::string config, rng;
-    this->build_filenames(traj, Params, config, smr, rng );
+    this->build_filenames(traj, Params, config, rng);
     this->check_filename(rng);
     this->check_filename(config);
 

Grid/qcd/hmc/checkpointers/ScidacCheckpointer.h

@@ -70,37 +70,19 @@ class ScidacHmcCheckpointer : public BaseHmcCheckpointer<Implementation> {
     }
   }
 
-  void TrajectoryComplete(int traj, 
+  void TrajectoryComplete(int traj, Field &U, GridSerialRNG &sRNG,
-			  ConfigurationBase<Field> &SmartConfig,
-			  GridSerialRNG &sRNG,
                           GridParallelRNG &pRNG) {
     if ((traj % Params.saveInterval) == 0) {
-      std::string config, rng,smr;
+      std::string config, rng;
-      this->build_filenames(traj, Params, config, smr, rng);
+      this->build_filenames(traj, Params, config, rng);
-      GridBase *grid = SmartConfig.get_U(false).Grid();
+      GridBase *grid = U.Grid();
       uint32_t nersc_csum,scidac_csuma,scidac_csumb;
       BinaryIO::writeRNG(sRNG, pRNG, rng, 0,nersc_csum,scidac_csuma,scidac_csumb);
-      std::cout << GridLogMessage << "Written Binary RNG " << rng
+      ScidacWriter _ScidacWriter(grid->IsBoss());
-                << " checksum " << std::hex 
+      _ScidacWriter.open(config);
-		<< nersc_csum   <<"/"
+      _ScidacWriter.writeScidacFieldRecord(U, MData);
-		<< scidac_csuma   <<"/"
+      _ScidacWriter.close();
-		<< scidac_csumb 
-		<< std::dec << std::endl;
 
-
-      {
-	ScidacWriter _ScidacWriter(grid->IsBoss());
-	_ScidacWriter.open(config);
-	_ScidacWriter.writeScidacFieldRecord(SmartConfig.get_U(false), MData);
-	_ScidacWriter.close();
-      }
-      
-      if ( Params.saveSmeared ) {
-	ScidacWriter _ScidacWriter(grid->IsBoss());
-	_ScidacWriter.open(smr);
-	_ScidacWriter.writeScidacFieldRecord(SmartConfig.get_U(true), MData);
-	_ScidacWriter.close();
-      }
       std::cout << GridLogMessage << "Written Scidac Configuration on " << config << std::endl;
     }
   };

Grid/qcd/hmc/integrators/Integrator.h

@@ -66,7 +66,6 @@ public:
 template <class FieldImplementation_, class SmearingPolicy, class RepresentationPolicy>
 class Integrator {
 protected:
-public:
   typedef FieldImplementation_ FieldImplementation;
   typedef typename FieldImplementation::Field MomentaField;  //for readability
   typedef typename FieldImplementation::Field Field;
@@ -87,8 +86,6 @@ public:
 
   const ActionSet<Field, RepresentationPolicy> as;
 
-  ActionSet<Field,RepresentationPolicy> LevelForces;
-  
   //Get a pointer to a shared static instance of the "do-nothing" momentum filter to serve as a default
   static MomentumFilterBase<MomentaField> const* getDefaultMomFilter(){ 
     static MomentumFilterNone<MomentaField> filter;
@@ -99,6 +96,7 @@ public:
   {
     t_P[level] += ep;
     update_P(P, U, level, ep);
+
     std::cout << GridLogIntegrator << "[" << level << "] P " << " dt " << ep << " : t_P " << t_P[level] << std::endl;
   }
 
@@ -126,33 +124,35 @@ public:
     // input U actually not used in the fundamental case
     // Fundamental updates, include smearing
 
-    assert(as.size()==LevelForces.size());
-    
-    Field level_force(U.Grid()); level_force =Zero();
     for (int a = 0; a < as[level].actions.size(); ++a) {
 
       double start_full = usecond();
       Field force(U.Grid());
       conformable(U.Grid(), Mom.Grid());
 
+      Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
       double start_force = usecond();
 
+      std::cout << GridLogMessage << "AuditForce["<<level<<"]["<<a<<"] before"<<std::endl;
+      
       as[level].actions.at(a)->deriv_timer_start();
-      as[level].actions.at(a)->deriv(Smearer, force);  // deriv should NOT include Ta
+      as[level].actions.at(a)->deriv(Us, force);  // deriv should NOT include Ta
       as[level].actions.at(a)->deriv_timer_stop();
 
+      std::cout << GridLogMessage << "AuditForce["<<level<<"]["<<a<<"] after"<<std::endl;
+
+      std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
       auto name = as[level].actions.at(a)->action_name();
+      if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
 
       force = FieldImplementation::projectForce(force); // Ta for gauge fields
       double end_force = usecond();
-      
+
+      //      DumpSliceNorm("force ",force,Nd-1);
       MomFilter->applyFilter(force);
-
       std::cout << GridLogIntegrator << " update_P : Level [" << level <<"]["<<a <<"] "<<name<<" dt "<<ep<<  std::endl;
-
+      DumpSliceNorm("force filtered ",force,Nd-1);
-      // track the total
+      
-      level_force = level_force+force;
-
       Real force_abs   = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm.  nb. norm2(latt) = \sum_x norm2(latt[x]) 
       Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR;    
 
@@ -175,16 +175,6 @@ public:
 
     }
 
-    {
-      // total force
-      Real force_abs   = std::sqrt(norm2(level_force)/U.Grid()->gSites()); //average per-site norm.  nb. norm2(latt) = \sum_x norm2(latt[x]) 
-      Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR;    
-
-      Real force_max   = std::sqrt(maxLocalNorm2(level_force));
-      Real impulse_max = force_max * ep * HMC_MOMENTUM_DENOMINATOR;    
-      LevelForces[level].actions.at(0)->deriv_log(force_abs,force_max,impulse_abs,impulse_max);
-    }
-
     // Force from the other representations
     as[level].apply(update_P_hireps, Representations, Mom, U, ep);
 
@@ -234,16 +224,6 @@ public:
 
     //Default the momentum filter to "do-nothing"
     MomFilter = getDefaultMomFilter();
-
-    for (int level = 0; level < as.size(); ++level) {
-      int multiplier = as.at(level).multiplier;
-      ActionLevel<Field, RepresentationPolicy> * Level = new ActionLevel<Field, RepresentationPolicy>(multiplier);
-      Level->push_back(new EmptyAction<Field>); 
-      LevelForces.push_back(*Level);
-      // does it copy by value or reference??
-      // - answer it copies by value, BUT the action level contains a reference that is NOT updated.
-      // Unsafe code in Guido's area
-    }
   };
 
   virtual ~Integrator() {}
@@ -261,14 +241,10 @@ public:
 
   void reset_timer(void)
   {
-    assert(as.size()==LevelForces.size());
     for (int level = 0; level < as.size(); ++level) {
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
         as[level].actions.at(actionID)->reset_timer();
       }
-      int actionID=0;
-      assert(LevelForces.at(level).actions.size()==1);
-      LevelForces.at(level).actions.at(actionID)->reset_timer();
     }
   }
   void print_timer(void)
@@ -330,16 +306,6 @@ public:
 		  <<" calls "     << as[level].actions.at(actionID)->deriv_num
 		  << std::endl;
       }
-      int actionID=0;
-      std::cout << GridLogMessage 
-		  << LevelForces[level].actions.at(actionID)->action_name()
-		  <<"["<<level<<"]["<< actionID<<"] :\n\t\t "
-		  <<" force max " << LevelForces[level].actions.at(actionID)->deriv_max_average()
-		  <<" norm "      << LevelForces[level].actions.at(actionID)->deriv_norm_average()
-		  <<" Fdt max  "  << LevelForces[level].actions.at(actionID)->Fdt_max_average()
-		  <<" Fdt norm "  << LevelForces[level].actions.at(actionID)->Fdt_norm_average()
-		  <<" calls "     << LevelForces[level].actions.at(actionID)->deriv_num
-		  << std::endl;
     }
     std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
   }
@@ -361,13 +327,6 @@ public:
 	std::cout << as[level].actions.at(actionID)->LogParameters();
       }
     }
-    std::cout << " [Integrator] Total Force loggers: "<< LevelForces.size() <<std::endl;
-    for (int level = 0; level < LevelForces.size(); ++level) {
-      std::cout << GridLogMessage << "[Integrator] ---- Level: "<< level << std::endl;
-      for (int actionID = 0; actionID < LevelForces[level].actions.size(); ++actionID) {
-	std::cout << GridLogMessage << "["<< LevelForces[level].actions.at(actionID)->action_name() << "] ID: " << actionID << std::endl;
-      }
-    }
     std::cout << GridLogMessage << ":::::::::::::::::::::::::::::::::::::::::"<< std::endl;
   }
 
@@ -418,9 +377,14 @@ public:
 	auto name = as[level].actions.at(actionID)->action_name();
         std::cout << GridLogMessage << "refresh [" << level << "][" << actionID << "] "<<name << std::endl;
 
+        Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+
+	std::cout << GridLogMessage << "AuditRefresh["<<level<<"]["<<actionID<<"] before"<<std::endl;
+
 	as[level].actions.at(actionID)->refresh_timer_start();
-        as[level].actions.at(actionID)->refresh(Smearer, sRNG, pRNG);
+        as[level].actions.at(actionID)->refresh(Us, sRNG, pRNG);
 	as[level].actions.at(actionID)->refresh_timer_stop();
+	std::cout << GridLogMessage << "AuditRefresh["<<level<<"]["<<actionID<<"] after"<<std::endl;
 
       }
 
@@ -449,7 +413,6 @@ public:
   RealD S(Field& U) 
   {  // here also U not used
 
-    assert(as.size()==LevelForces.size());
     std::cout << GridLogIntegrator << "Integrator action\n";
 
     RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom
@@ -462,9 +425,10 @@ public:
 
         // get gauge field from the SmearingPolicy and
         // based on the boolean is_smeared in actionID
+        Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] action eval " << std::endl;
 	        as[level].actions.at(actionID)->S_timer_start();
-        Hterm = as[level].actions.at(actionID)->S(Smearer);
+        Hterm = as[level].actions.at(actionID)->S(Us);
    	        as[level].actions.at(actionID)->S_timer_stop();
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] H = " << Hterm << std::endl;
         H += Hterm;
@@ -505,11 +469,12 @@ public:
       for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) {
         // get gauge field from the SmearingPolicy and
         // based on the boolean is_smeared in actionID
+        Field& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] action eval " << std::endl;
+	        as[level].actions.at(actionID)->S_timer_start();
 
-	as[level].actions.at(actionID)->S_timer_start();
+        Hterm = as[level].actions.at(actionID)->Sinitial(Us);
-        Hterm = as[level].actions.at(actionID)->S(Smearer);
+   	        as[level].actions.at(actionID)->S_timer_stop();
-	as[level].actions.at(actionID)->S_timer_stop();
 
         std::cout << GridLogMessage << "S [" << level << "][" << actionID << "] H = " << Hterm << std::endl;
         H += Hterm;

Grid/qcd/observables/hmc_observable.h

@@ -34,13 +34,6 @@ NAMESPACE_BEGIN(Grid);
 template <class Field>
 class HmcObservable {
  public:
-  virtual void TrajectoryComplete(int traj,
-                                  ConfigurationBase<Field> &SmartConfig,
-                                  GridSerialRNG &sRNG,
-                                  GridParallelRNG &pRNG)
-  {
-    TrajectoryComplete(traj,SmartConfig.get_U(false),sRNG,pRNG); // Unsmeared observable
-  };
   virtual void TrajectoryComplete(int traj,
                                   Field &U,
                                   GridSerialRNG &sRNG,

Grid/qcd/observables/plaquette.h

@@ -42,18 +42,6 @@ public:
   // necessary for HmcObservable compatibility
   typedef typename Impl::Field Field;
 
-  virtual void TrajectoryComplete(int traj,
-                                  ConfigurationBase<Field> &SmartConfig,
-                                  GridSerialRNG &sRNG,
-                                  GridParallelRNG &pRNG)
-  {
-    std::cout << GridLogMessage << "+++++++++++++++++++"<<std::endl;
-    std::cout << GridLogMessage << "Unsmeared plaquette"<<std::endl;
-    TrajectoryComplete(traj,SmartConfig.get_U(false),sRNG,pRNG); // Unsmeared observable
-    std::cout << GridLogMessage << "Smeared plaquette"<<std::endl;
-    TrajectoryComplete(traj,SmartConfig.get_U(true),sRNG,pRNG); // Unsmeared observable
-    std::cout << GridLogMessage << "+++++++++++++++++++"<<std::endl;
-  };
   void TrajectoryComplete(int traj,
                           Field &U,
                           GridSerialRNG &sRNG,

Grid/qcd/representations/fundamental.h

@@ -13,7 +13,7 @@ NAMESPACE_BEGIN(Grid);
  * Empty since HMC updates already the fundamental representation 
  */
 
-template <int ncolour, class group_name>
+template <int ncolour>
 class FundamentalRep {
 public:
   static const int Dimension = ncolour;
@@ -21,7 +21,7 @@ public:
 
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields
-  typedef typename GaugeGroup<ncolour,group_name>::LatticeMatrix LatticeMatrix;
+  typedef typename SU<ncolour>::LatticeMatrix LatticeMatrix;
   typedef LatticeGaugeField LatticeField;
   
   explicit FundamentalRep(GridBase* grid) {} //do nothing
@@ -45,8 +45,7 @@ public:
     
 
   
-typedef	 FundamentalRep<Nc,GroupName::SU> FundamentalRepresentation;
+typedef	 FundamentalRep<Nc> FundamentalRepresentation;
-typedef	 FundamentalRep<Nc,GroupName::Sp> SpFundamentalRepresentation;
 
 NAMESPACE_END(Grid);  
 

Grid/qcd/representations/two_index.h

@@ -20,14 +20,14 @@ NAMESPACE_BEGIN(Grid);
  * in the SUnTwoIndex.h file
  */
 
-template <int ncolour, TwoIndexSymmetry S, class group_name = GroupName::SU>
+template <int ncolour, TwoIndexSymmetry S>
 class TwoIndexRep {
 public:
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields
-  typedef typename GaugeGroupTwoIndex<ncolour, S, group_name>::LatticeTwoIndexMatrix LatticeMatrix;
+  typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexMatrix LatticeMatrix;
-  typedef typename GaugeGroupTwoIndex<ncolour, S, group_name>::LatticeTwoIndexField LatticeField;
+  typedef typename SU_TwoIndex<ncolour, S>::LatticeTwoIndexField LatticeField;
-  static const int Dimension = GaugeGroupTwoIndex<ncolour,S,group_name>::Dimension;
+  static const int Dimension = ncolour * (ncolour + S) / 2;
   static const bool isFundamental = false;
 
   LatticeField U;
@@ -43,10 +43,10 @@ public:
     U = Zero();
     LatticeColourMatrix tmp(Uin.Grid());
 
-    Vector<typename GaugeGroup<ncolour,group_name>::Matrix> eij(Dimension);
+    Vector<typename SU<ncolour>::Matrix> eij(Dimension);
 
     for (int a = 0; a < Dimension; a++)
-      GaugeGroupTwoIndex<ncolour, S, group_name>::base(a, eij[a]);
+      SU_TwoIndex<ncolour, S>::base(a, eij[a]);
 
     for (int mu = 0; mu < Nd; mu++) {
       auto Uin_mu = peekLorentz(Uin, mu);
@@ -71,7 +71,7 @@ public:
 
       out_mu = Zero();
 
-      typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector h(in.Grid());
+      typename SU<ncolour>::LatticeAlgebraVector h(in.Grid());
       projectOnAlgebra(h, in_mu, double(Nc + 2 * S));  // factor T(r)/T(fund)
       FundamentalLieAlgebraMatrix(h, out_mu);          // apply scale only once
       pokeLorentz(out, out_mu, mu);
@@ -80,23 +80,20 @@ public:
   }
 
 private:
-  void projectOnAlgebra(typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h_out,
+  void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out,
                         const LatticeMatrix &in, Real scale = 1.0) const {
-    GaugeGroupTwoIndex<ncolour, S,group_name>::projectOnAlgebra(h_out, in, scale);
+    SU_TwoIndex<ncolour, S>::projectOnAlgebra(h_out, in, scale);
   }
 
   void FundamentalLieAlgebraMatrix(
-				   typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h,
+				   typename SU<ncolour>::LatticeAlgebraVector &h,
-				   typename GaugeGroup<ncolour, group_name>::LatticeMatrix &out, Real scale = 1.0) const {
+				   typename SU<ncolour>::LatticeMatrix &out, Real scale = 1.0) const {
-    GaugeGroup<ncolour,group_name>::FundamentalLieAlgebraMatrix(h, out, scale);
+    SU<ncolour>::FundamentalLieAlgebraMatrix(h, out, scale);
   }
 };
 
-typedef TwoIndexRep<Nc, Symmetric, GroupName::SU> TwoIndexSymmetricRepresentation;
+typedef TwoIndexRep<Nc, Symmetric> TwoIndexSymmetricRepresentation;
-typedef TwoIndexRep<Nc, AntiSymmetric, GroupName::SU> TwoIndexAntiSymmetricRepresentation;
+typedef TwoIndexRep<Nc, AntiSymmetric> TwoIndexAntiSymmetricRepresentation;
-
-typedef TwoIndexRep<Nc, Symmetric, GroupName::Sp> SpTwoIndexSymmetricRepresentation;
-typedef TwoIndexRep<Nc, AntiSymmetric, GroupName::Sp> SpTwoIndexAntiSymmetricRepresentation;
 
 NAMESPACE_END(Grid);
 

Grid/qcd/smearing/GaugeConfiguration.h

@@ -7,27 +7,26 @@
 
 NAMESPACE_BEGIN(Grid);
 
-
 //trivial class for no smearing
 template< class Impl >
-class NoSmearing : public ConfigurationBase<typename Impl::Field>
+class NoSmearing
 {
 public:
   INHERIT_FIELD_TYPES(Impl);
 
-  Field* ThinLinks;
+  Field* ThinField;
 
-  NoSmearing(): ThinLinks(NULL) {}
+  NoSmearing(): ThinField(NULL) {}
 
-  virtual void set_Field(Field& U) { ThinLinks = &U; }
+  void set_Field(Field& U) { ThinField = &U; }
 
-  virtual void smeared_force(Field&) {}
+  void smeared_force(Field&) const {}
 
-  virtual Field& get_SmearedU() { return *ThinLinks; }
+  Field& get_SmearedU() { return *ThinField; }
 
-  virtual Field &get_U(bool smeared = false)
+  Field &get_U(bool smeared = false)
   {
-    return *ThinLinks;
+    return *ThinField;
   }
 };
 
@@ -43,24 +42,19 @@ public:
   It stores a list of smeared configurations.
 */
 template <class Gimpl>
-class SmearedConfiguration : public ConfigurationBase<typename Gimpl::Field>
+class SmearedConfiguration
 {
 public:
   INHERIT_GIMPL_TYPES(Gimpl);
 
-protected:
+private:
   const unsigned int smearingLevels;
   Smear_Stout<Gimpl> *StoutSmearing;
   std::vector<GaugeField> SmearedSet;
-public:
+
-  GaugeField*  ThinLinks; /* Pointer to the thin links configuration */ // move to base???
-protected:
-  
   // Member functions
   //====================================================================
-
+  void fill_smearedSet(GaugeField &U)
-  // Overridden in masked version
-  virtual void fill_smearedSet(GaugeField &U)
   {
     ThinLinks = &U;  // attach the smearing routine to the field U
 
@@ -88,10 +82,9 @@ protected:
       }
     }
   }
-
+  //====================================================================
-  //overridden in masked verson
+  GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
-  virtual GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime,
+                                  const GaugeField& GaugeK) const 
-					  const GaugeField& GaugeK) const 
   {
     GridBase* grid = GaugeK.Grid();
     GaugeField C(grid), SigmaK(grid), iLambda(grid);
@@ -220,6 +213,8 @@ protected:
 
   //====================================================================
 public:
+  GaugeField*
+      ThinLinks; /* Pointer to the thin links configuration */
 
   /* Standard constructor */
   SmearedConfiguration(GridCartesian* UGrid, unsigned int Nsmear,
@@ -235,7 +230,7 @@ public:
     : smearingLevels(0), StoutSmearing(nullptr), SmearedSet(), ThinLinks(NULL) {}
 
   // attach the smeared routines to the thin links U and fill the smeared set
-  virtual void set_Field(GaugeField &U)
+  void set_Field(GaugeField &U)
   {
     double start = usecond();
     fill_smearedSet(U);
@@ -245,7 +240,7 @@ public:
   }
 
   //====================================================================
-  virtual void smeared_force(GaugeField &SigmaTilde) 
+  void smeared_force(GaugeField &SigmaTilde) const
   {
     if (smearingLevels > 0)
     {
@@ -272,16 +267,14 @@ public:
       }
       double end = usecond();
       double time = (end - start)/ 1e3;
-      std::cout << GridLogMessage << " GaugeConfiguration: Smeared Force chain rule took " << time << " ms" << std::endl;
+      std::cout << GridLogMessage << "Smearing force in " << time << " ms" << std::endl;  
     }  // if smearingLevels = 0 do nothing
-    SigmaTilde=Gimpl::projectForce(SigmaTilde); // Ta
-      
   }
   //====================================================================
 
-  virtual GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; }
+  GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; }
 
-  virtual GaugeField &get_U(bool smeared = false)
+  GaugeField &get_U(bool smeared = false)
   {
     // get the config, thin links by default
     if (smeared)

Grid/qcd/smearing/HISQSmearing.h

@@ -1,389 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/smearing/HISQSmearing.h
-
-Copyright (C) 2023
-
-Author: D. A. Clarke <clarke.davida@gmail.com> 
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-/*
-    @file HISQSmearing.h
-    @brief Declares classes related to HISQ smearing 
-*/
-
-
-#pragma once
-#include <Grid/Grid.h>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-
-NAMESPACE_BEGIN(Grid);
-
-
-// TODO: find a way to fold this into the stencil header. need to access grid to get
-// Nd, since you don't want to inherit from QCD.h
-/*!  @brief append arbitrary shift path to shifts */
-template<typename... Args>
-void appendShift(std::vector<Coordinate>& shifts, int dir, Args... args) {
-    Coordinate shift(Nd,0);
-    generalShift(shift, dir, args...); 
-    // push_back creates an element at the end of shifts and
-    // assigns the data in the argument to it.
-    shifts.push_back(shift);
-}
-
-
-/*!  @brief figure out the stencil index from mu and nu */
-accelerator_inline int stencilIndex(int mu, int nu) {
-    // Nshifts depends on how you built the stencil
-    int Nshifts = 6;
-    return Nshifts*nu + Nd*Nshifts*mu;
-}
-
-
-/*!  @brief structure holding the link treatment */
-struct SmearingParameters{
-    SmearingParameters(){}
-    Real c_1;               // 1 link
-    Real c_naik;            // Naik term
-    Real c_3;               // 3 link
-    Real c_5;               // 5 link
-    Real c_7;               // 7 link
-    Real c_lp;              // 5 link Lepage
-    SmearingParameters(Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) 
-        : c_1(c1),
-          c_naik(cnaik),
-          c_3(c3),
-          c_5(c5),
-          c_7(c7),
-          c_lp(clp){}
-};
-
-
-/*!  @brief create fat links from link variables */
-template<class Gimpl> 
-class Smear_HISQ : public Gimpl {
-
-private:
-    GridCartesian* const _grid;
-    SmearingParameters _linkTreatment;
-
-public:
-
-    INHERIT_GIMPL_TYPES(Gimpl);
-    typedef typename Gimpl::GaugeField     GF;
-    typedef typename Gimpl::GaugeLinkField LF;
-    typedef typename Gimpl::ComplexField   CF;
-
-    // Don't allow default values here.
-    Smear_HISQ(GridCartesian* grid, Real c1, Real cnaik, Real c3, Real c5, Real c7, Real clp) 
-        : _grid(grid), 
-          _linkTreatment(c1,cnaik,c3,c5,c7,clp) {
-        assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3");
-        assert(Nd == 4 && "HISQ smearing only defined for Nd==4");
-    }
-
-    // Allow to pass a pointer to a C-style, double array for MILC convenience
-    Smear_HISQ(GridCartesian* grid, double* coeff) 
-        : _grid(grid), 
-          _linkTreatment(coeff[0],coeff[1],coeff[2],coeff[3],coeff[4],coeff[5]) {
-        assert(Nc == 3 && "HISQ smearing currently implemented only for Nc==3");
-        assert(Nd == 4 && "HISQ smearing only defined for Nd==4");
-    }
-
-    ~Smear_HISQ() {}
-
-    // Intent: OUT--u_smr, u_naik
-    //          IN--u_thin
-    void smear(GF& u_smr, GF& u_naik, GF& u_thin) const {
-
-        SmearingParameters lt = this->_linkTreatment;
-        auto grid = this->_grid;
-
-        // Create a padded cell of extra padding depth=1 and fill the padding.
-        int depth = 1;
-        PaddedCell Ghost(depth,grid);
-        GF Ughost = Ghost.Exchange(u_thin);
-
-        // This is where auxiliary N-link fields and the final smear will be stored. 
-        GF Ughost_fat(Ughost.Grid());
-        GF Ughost_3link(Ughost.Grid());
-        GF Ughost_5linkA(Ughost.Grid());
-        GF Ughost_5linkB(Ughost.Grid());
-
-        // mu-nu plane stencil. We allow mu==nu to make indexing the stencil easier,
-        // but these entries will not be used. 
-        std::vector<Coordinate> shifts;
-        for(int mu=0;mu<Nd;mu++)
-        for(int nu=0;nu<Nd;nu++) {
-            appendShift(shifts,mu);
-            appendShift(shifts,nu);
-            appendShift(shifts,shiftSignal::NO_SHIFT);
-            appendShift(shifts,mu,Back(nu));
-            appendShift(shifts,Back(nu));
-            appendShift(shifts,Back(mu));
-        }
-
-        // A GeneralLocalStencil has two indices: a site and stencil index 
-        GeneralLocalStencil gStencil(Ughost.Grid(),shifts);
-
-        // This is where contributions from the smearing get added together
-        Ughost_fat=Zero();
-
-        // This loop handles 3-, 5-, and 7-link constructs, minus Lepage and Naik.
-        for(int mu=0;mu<Nd;mu++) {
-
-            // TODO: This approach is slightly memory inefficient. It uses 25% extra memory 
-            Ughost_3link =Zero();
-            Ughost_5linkA=Zero();
-            Ughost_5linkB=Zero();
-
-            // Create the accessors
-            autoView(U_v       , Ughost       , AcceleratorRead);
-            autoView(U_fat_v   , Ughost_fat   , AcceleratorWrite);
-            autoView(U_3link_v , Ughost_3link , AcceleratorWrite);
-            autoView(U_5linkA_v, Ughost_5linkA, AcceleratorWrite);
-            autoView(U_5linkB_v, Ughost_5linkB, AcceleratorWrite);
-
-            // We infer some types that will be needed in the calculation.
-            typedef decltype(gStencil.GetEntry(0,0)) stencilElement;
-            typedef decltype(coalescedReadGeneralPermute(U_v[0](0),gStencil.GetEntry(0,0)->_permute,Nd)) U3matrix;
-
-            int Nsites = U_v.size();
-            auto gStencil_v = gStencil.View(); 
-
-            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 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++) {
-                    if(nu==mu) continue;
-                    int s = stencilIndex(mu,nu);
-
-                    // The stencil gives us support points in the mu-nu plane that we will use to
-                    // grab the links we need.
-                    SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
-                    SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
-                    SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
-                    SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
-                    SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
-                    SE5 = gStencil_v.GetEntry(s+5,site); int x_m_mu      = SE5->_offset;
-
-                    // When you're deciding whether to take an adjoint, the question is: how is the
-                    // stored link oriented compared to the one you want? If I imagine myself travelling
-                    // with the to-be-updated link, I have two possible, alternative 3-link paths I can
-                    // take, one starting by going to the left, the other starting by going to the right.
-                    U0 = coalescedReadGeneralPermute(U_v[x_p_mu     ](nu),SE0->_permute,Nd);
-                    U1 = coalescedReadGeneralPermute(U_v[x_p_nu     ](mu),SE1->_permute,Nd);
-                    U2 = coalescedReadGeneralPermute(U_v[x          ](nu),SE2->_permute,Nd);
-                    U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd);
-                    U4 = coalescedReadGeneralPermute(U_v[x_m_nu     ](mu),SE4->_permute,Nd);
-                    U5 = coalescedReadGeneralPermute(U_v[x_m_nu     ](nu),SE4->_permute,Nd);
-
-                    //  "left"          "right"
-                    W = U2*U1*adj(U0) + adj(U5)*U4*U3;
-
-                    // Save 3-link construct for later and add to smeared field.
-                    coalescedWrite(U_3link_v[x](nu), W);
-
-                    // The index operator (x) returns the coalesced read on GPU. The view [] index returns 
-                    // a reference to the vector object. The [x](mu) returns a reference to the densely 
-                    // packed (contiguous in memory) mu-th element of the vector object. On CPU, 
-                    // coalescedRead/Write is the identity mapping assigning vector object to vector object.
-                    // But on GPU it's non-trivial and maps scalar object to vector object and vice versa.
-                    coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_3*W);
-                }
-            })
-
-            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 5-link 
-                stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
-                U3matrix U0, U1, U2, U3, U4, U5, W;
-                int sigmaIndex = 0;
-                for(int nu=0;nu<Nd;nu++) {
-                    if(nu==mu) continue;
-                    int s = stencilIndex(mu,nu);
-                    for(int rho=0;rho<Nd;rho++) {
-                        if (rho == mu || rho == nu) continue;
-
-                        SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
-                        SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
-                        SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
-                        SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
-                        SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
-
-                        U0 = coalescedReadGeneralPermute(      U_v[x_p_mu     ](nu ),SE0->_permute,Nd);
-                        U1 = coalescedReadGeneralPermute(U_3link_v[x_p_nu     ](rho),SE1->_permute,Nd);
-                        U2 = coalescedReadGeneralPermute(      U_v[x          ](nu ),SE2->_permute,Nd);
-                        U3 = coalescedReadGeneralPermute(      U_v[x_p_mu_m_nu](nu ),SE3->_permute,Nd);
-                        U4 = coalescedReadGeneralPermute(U_3link_v[x_m_nu     ](rho),SE4->_permute,Nd);
-                        U5 = coalescedReadGeneralPermute(      U_v[x_m_nu     ](nu ),SE4->_permute,Nd);
-
-                        W  = U2*U1*adj(U0) + adj(U5)*U4*U3;
-
-                        if(sigmaIndex<3) {
-                            coalescedWrite(U_5linkA_v[x](rho), W);
-                        } else {
-                            coalescedWrite(U_5linkB_v[x](rho), W);
-                        }    
-
-                        coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_5*W);
-                        sigmaIndex++;
-                    }
-                }
-            })
-
-            accelerator_for(site,Nsites,Simd::Nsimd(),{ // ----------- 7-link
-                stencilElement SE0, SE1, SE2, SE3, SE4, SE5;
-                U3matrix U0, U1, U2, U3, U4, U5, W;
-                int sigmaIndex = 0;
-                for(int nu=0;nu<Nd;nu++) {
-                    if(nu==mu) continue;
-                    int s = stencilIndex(mu,nu);
-                    for(int rho=0;rho<Nd;rho++) {
-                        if (rho == mu || rho == nu) continue;
-
-                        SE0 = gStencil_v.GetEntry(s+0,site); int x_p_mu      = SE0->_offset;
-                        SE1 = gStencil_v.GetEntry(s+1,site); int x_p_nu      = SE1->_offset;
-                        SE2 = gStencil_v.GetEntry(s+2,site); int x           = SE2->_offset;
-                        SE3 = gStencil_v.GetEntry(s+3,site); int x_p_mu_m_nu = SE3->_offset;
-                        SE4 = gStencil_v.GetEntry(s+4,site); int x_m_nu      = SE4->_offset;
-
-                        U0 = coalescedReadGeneralPermute(U_v[x_p_mu](nu),SE0->_permute,Nd);
-                        if(sigmaIndex<3) {
-                            U1 = coalescedReadGeneralPermute(U_5linkB_v[x_p_nu](rho),SE1->_permute,Nd);
-                        } else {
-                            U1 = coalescedReadGeneralPermute(U_5linkA_v[x_p_nu](rho),SE1->_permute,Nd);
-                        }  
-                        U2 = coalescedReadGeneralPermute(U_v[x](nu),SE2->_permute,Nd);
-                        U3 = coalescedReadGeneralPermute(U_v[x_p_mu_m_nu](nu),SE3->_permute,Nd);
-                        if(sigmaIndex<3) {
-                            U4 = coalescedReadGeneralPermute(U_5linkB_v[x_m_nu](rho),SE4->_permute,Nd);
-                        } else {
-                            U4 = coalescedReadGeneralPermute(U_5linkA_v[x_m_nu](rho),SE4->_permute,Nd);
-                        }  
-                        U5 = coalescedReadGeneralPermute(U_v[x_m_nu](nu),SE4->_permute,Nd);
-
-                        W  = U2*U1*adj(U0) + adj(U5)*U4*U3;
-
-                        coalescedWrite(U_fat_v[x](mu), U_fat_v(x)(mu) + lt.c_7*W);
-                        sigmaIndex++;
-                    }
-                }
-            })
-
-        } // end mu loop
-
-        // c1, c3, c5, c7 construct contributions
-        u_smr = Ghost.Extract(Ughost_fat) + lt.c_1*u_thin;
-
-        // Load up U and V std::vectors to access thin and smeared links.
-        std::vector<LF> U(Nd, grid);
-        std::vector<LF> V(Nd, grid);
-        std::vector<LF> Vnaik(Nd, grid);
-        for (int mu = 0; mu < Nd; mu++) {
-            U[mu] = PeekIndex<LorentzIndex>(u_thin, mu);
-            V[mu] = PeekIndex<LorentzIndex>(u_smr, mu);
-        }
-
-        for(int mu=0;mu<Nd;mu++) {
-
-            // Naik
-            Vnaik[mu] = lt.c_naik*Gimpl::CovShiftForward(U[mu],mu,
-                                    Gimpl::CovShiftForward(U[mu],mu,
-                                      Gimpl::CovShiftIdentityForward(U[mu],mu)));
-
-            // LePage
-            for (int nu_h=1;nu_h<Nd;nu_h++) {
-                int nu=(mu+nu_h)%Nd;
-                                // nu, nu, mu, Back(nu), Back(nu)
-                V[mu] = V[mu] + lt.c_lp*Gimpl::CovShiftForward(U[nu],nu,
-                                          Gimpl::CovShiftForward(U[nu],nu,
-                                            Gimpl::CovShiftForward(U[mu],mu,
-                                              Gimpl::CovShiftBackward(U[nu],nu,
-                                                Gimpl::CovShiftIdentityBackward(U[nu],nu)))))
-                                // Back(nu), Back(nu), mu, nu, nu
-                              + lt.c_lp*Gimpl::CovShiftBackward(U[nu],nu,
-                                          Gimpl::CovShiftBackward(U[nu],nu,
-                                            Gimpl::CovShiftForward(U[mu],mu,
-                                              Gimpl::CovShiftForward(U[nu],nu,
-                                                Gimpl::CovShiftIdentityForward(U[nu],nu)))));
-            }
-        }
-
-        // Put V back into u_smr.
-        for (int mu = 0; mu < Nd; mu++) {
-            PokeIndex<LorentzIndex>(u_smr , V[mu]    , mu);
-            PokeIndex<LorentzIndex>(u_naik, Vnaik[mu], mu);
-        }
-    };
-
-
-    // Intent: OUT--u_proj
-    //          IN--u_mu
-    void projectU3(GF& u_proj, GF& u_mu) const {
-
-        auto grid = this->_grid;
-
-        LF V(grid), Q(grid), sqrtQinv(grid), id_3(grid), diff(grid);
-        CF c0(grid), c1(grid), c2(grid), g0(grid), g1(grid), g2(grid), S(grid), R(grid), theta(grid), 
-           u(grid), v(grid), w(grid), den(grid), f0(grid), f1(grid), f2(grid);
-
-        // Follow MILC 10.1103/PhysRevD.82.074501, eqs (B2-B3) and (C1-C8)
-        for (int mu = 0; mu < Nd; mu++) {
-            V  = PeekIndex<LorentzIndex>(u_mu, mu);
-            Q  = adj(V)*V;
-            c0 =        real(trace(Q));
-            c1 = (1/2.)*real(trace(Q*Q));
-            c2 = (1/3.)*real(trace(Q*Q*Q));
-            S  = (1/3.)*c1-(1/18.)*c0*c0;
-            if (norm2(S)<1e-28) {
-                g0 = (1/3.)*c0; g1 = g0; g2 = g1;
-            } else {
-                R     = (1/2.)*c2-(1/3. )*c0*c1+(1/27.)*c0*c0*c0;
-                theta = acos(R*pow(S,-1.5));
-                g0    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta-2*M_PI/3.);
-                g1    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta          );
-                g2    = (1/3.)*c0+2.*sqrt(S)*cos((1/3.)*theta+2*M_PI/3.);
-            }
-//            if (fabs(Q.determinant()/(g0*g1*g2)-1.0) > 1e-5) { SVD }
-            u     = sqrt(g0) + sqrt(g1) + sqrt(g2);
-            v     = sqrt(g0*g1) + sqrt(g0*g2) + sqrt(g1*g2);
-            w     = sqrt(g0*g1*g2);
-            den   = w*(u*v-w);
-            f0    = (-w*(u*u+v)+u*v*v)/den;
-            f1    = (-w-u*u*u+2.*u*v)/den;
-            f2    = u/den;
-            id_3  = 1.;
-
-            sqrtQinv = f0*id_3 + f1*Q + f2*Q*Q;
-
-            PokeIndex<LorentzIndex>(u_proj, V*sqrtQinv, mu);
-        }
-    };
-
-
-//    void derivative(const GaugeField& Gauge) const {
-//    };
-};
-
-
-NAMESPACE_END(Grid);

Grid/qcd/smearing/JacobianAction.h

@@ -1,87 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/action/gauge/JacobianAction.h
-
-Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution
-directory
-*************************************************************************************/
-			   /*  END LEGAL */
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-////////////////////////////////////////////////////////////////////////
-// Jacobian Action .. 
-////////////////////////////////////////////////////////////////////////
-template <class Gimpl>
-class JacobianAction : public Action<typename Gimpl::GaugeField> {
-public:  
-  INHERIT_GIMPL_TYPES(Gimpl);
-
-  SmearedConfigurationMasked<Gimpl> * smearer;
-  /////////////////////////// constructors
-  explicit JacobianAction(SmearedConfigurationMasked<Gimpl> * _smearer ) { smearer=_smearer;};
-
-  virtual std::string action_name() {return "JacobianAction";}
-
-  virtual std::string LogParameters(){
-    std::stringstream sstream;
-    sstream << GridLogMessage << "[JacobianAction] " << std::endl;
-    return sstream.str();
-  }
-
-  //////////////////////////////////
-  // Usual cases are not used
-  //////////////////////////////////
-  virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){ assert(0);};
-  virtual RealD S(const GaugeField &U) { assert(0); }
-  virtual void deriv(const GaugeField &U, GaugeField &dSdU) { assert(0);  }
-
-  //////////////////////////////////
-  // Functions of smart configs only
-  //////////////////////////////////
-  virtual void refresh(ConfigurationBase<GaugeField> & U, GridSerialRNG &sRNG, GridParallelRNG& pRNG)
-  {
-    return;
-  }
-  virtual RealD S(ConfigurationBase<GaugeField>& U)
-  {
-    // det M = e^{ - ( - logDetM) }
-    assert( &U == smearer );
-    return -smearer->logDetJacobian();
-  }
-  virtual RealD Sinitial(ConfigurationBase<GaugeField>& U) 
-  {
-    return S(U);
-  }
-  virtual void deriv(ConfigurationBase<GaugeField>& U, GaugeField& dSdU)
-  {
-    assert( &U == smearer );
-    smearer->logDetJacobianForce(dSdU);
-  }
-
-private:
- };
-
-NAMESPACE_END(Grid);
-

Grid/qcd/smearing/Smearing.h

@@ -5,5 +5,4 @@
 #include <Grid/qcd/smearing/StoutSmearing.h>
 #include <Grid/qcd/smearing/GaugeConfiguration.h>
 #include <Grid/qcd/smearing/WilsonFlow.h>
-#include <Grid/qcd/smearing/HISQSmearing.h>
 

Grid/qcd/smearing/StoutSmearing.h

@@ -40,9 +40,7 @@ template <class Gimpl>
 class Smear_Stout : public Smear<Gimpl> {
  private:
   int OrthogDim = -1;
-public:
   const std::vector<double> SmearRho;
-private:
   // Smear<Gimpl>* ownership semantics:
   //    Smear<Gimpl>* passed in to constructor are owned by caller, so we don't delete them here
   //    Smear<Gimpl>* created within constructor need to be deleted as part of the destructor
@@ -69,7 +67,7 @@ public:
   /*! Construct stout smearing object from explicitly specified rho matrix */
   Smear_Stout(const std::vector<double>& rho_)
     : OwnedBase{new Smear_APE<Gimpl>(rho_)}, SmearBase{OwnedBase.get()} {
-    std::cout << GridLogDebug << "Stout smearing constructor : Smear_Stout(const std::vector<double>& " << rho_ << " )" << std::endl;
+    std::cout << GridLogDebug << "Stout smearing constructor : Smear_Stout(const std::vector<double>& " << rho_ << " )" << std::endl
     assert(Nc == 3 && "Stout smearing currently implemented only for Nc==3");
     }
 

Grid/qcd/utils/CovariantCshift.h

@@ -37,14 +37,13 @@ NAMESPACE_BEGIN(Grid);
 // Make these members of an Impl class for BC's.
 
 namespace PeriodicBC { 
-  //Out(x) = Link(x)*field(x+mu)
+
   template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
 									   int mu,
 									   const Lattice<covariant> &field)
   {
     return Link*Cshift(field,mu,1);// moves towards negative mu
   }
-  //Out(x) = Link^dag(x-mu)*field(x-mu)
   template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link, 
 									    int mu,
 									    const Lattice<covariant> &field)
@@ -53,19 +52,19 @@ namespace PeriodicBC {
     tmp = adj(Link)*field;
     return Cshift(tmp,mu,-1);// moves towards positive mu
   }
-  //Out(x) = Link^dag(x-mu)
+
   template<class gauge> Lattice<gauge>
   CovShiftIdentityBackward(const Lattice<gauge> &Link, int mu) 
   {
     return Cshift(adj(Link), mu, -1);
   }
-  //Out(x) = Link(x)
+
   template<class gauge> Lattice<gauge>
   CovShiftIdentityForward(const Lattice<gauge> &Link, int mu)
   {
     return Link;
   }
-  //Link(x) = Link(x+mu)
+
   template<class gauge> Lattice<gauge>
   ShiftStaple(const Lattice<gauge> &Link, int mu)
   {

Grid/qcd/utils/GaugeGroup.h

@@ -1,528 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid
-
-Source file: ./lib/qcd/utils/GaugeGroup.h
-
-Copyright (C) 2015
-
-Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-Author: neo <cossu@post.kek.jp>
-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 */
-#ifndef QCD_UTIL_GAUGEGROUP_H
-#define QCD_UTIL_GAUGEGROUP_H
-
-// Important detail: nvcc requires all template parameters to have names.
-// This is the only reason why the second template parameter has a name.
-#define ONLY_IF_SU                                                       \
-  typename dummy_name = group_name,                                      \
-           typename named_dummy = std::enable_if_t <                                 \
-                          std::is_same<dummy_name, group_name>::value && \
-                      is_su<dummy_name>::value >
-
-#define ONLY_IF_Sp                                                       \
-  typename dummy_name = group_name,                                      \
-           typename named_dummy = std::enable_if_t <                                 \
-                          std::is_same<dummy_name, group_name>::value && \
-                      is_sp<dummy_name>::value >
-
-NAMESPACE_BEGIN(Grid);
-namespace GroupName {
-class SU {};
-class Sp {};
-}  // namespace GroupName
-
-template <typename group_name>
-struct is_su {
-  static const bool value = false;
-};
-
-template <>
-struct is_su<GroupName::SU> {
-  static const bool value = true;
-};
-
-template <typename group_name>
-struct is_sp {
-  static const bool value = false;
-};
-
-template <>
-struct is_sp<GroupName::Sp> {
-  static const bool value = true;
-};
-
-template <typename group_name>
-constexpr int compute_adjoint_dimension(int ncolour);
-
-template <>
-constexpr int compute_adjoint_dimension<GroupName::SU>(int ncolour) {
-  return ncolour * ncolour - 1;
-}
-
-template <>
-constexpr int compute_adjoint_dimension<GroupName::Sp>(int ncolour) {
-  return ncolour / 2 * (ncolour + 1);
-}
-
-template <int ncolour, class group_name>
-class GaugeGroup {
- public:
-  static const int Dimension = ncolour;
-  static const int AdjointDimension =
-      compute_adjoint_dimension<group_name>(ncolour);
-  static const int AlgebraDimension =
-      compute_adjoint_dimension<group_name>(ncolour);
-
-  template <typename vtype>
-  using iSU2Matrix = iScalar<iScalar<iMatrix<vtype, 2> > >;
-  template <typename vtype>
-  using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
-  template <typename vtype>
-  using iAlgebraVector = iScalar<iScalar<iVector<vtype, AdjointDimension> > >;
-  template <typename vtype>
-  using iSUnAlgebraMatrix =
-    iScalar<iScalar<iMatrix<vtype, AdjointDimension> > >;
-  static int su2subgroups(void) { return su2subgroups(group_name()); }
-
-  //////////////////////////////////////////////////////////////////////////////////////////////////
-  // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix,
-  // SU<2>::LatticeMatrix etc...
-  //////////////////////////////////////////////////////////////////////////////////////////////////
-  typedef iGroupMatrix<Complex> Matrix;
-  typedef iGroupMatrix<ComplexF> MatrixF;
-  typedef iGroupMatrix<ComplexD> MatrixD;
-
-  typedef iGroupMatrix<vComplex> vMatrix;
-  typedef iGroupMatrix<vComplexF> vMatrixF;
-  typedef iGroupMatrix<vComplexD> vMatrixD;
-
-  // For the projectors to the algebra
-  // these should be real...
-  // keeping complex for consistency with the SIMD vector types
-  typedef iAlgebraVector<Complex> AlgebraVector;
-  typedef iAlgebraVector<ComplexF> AlgebraVectorF;
-  typedef iAlgebraVector<ComplexD> AlgebraVectorD;
-
-  typedef iAlgebraVector<vComplex> vAlgebraVector;
-  typedef iAlgebraVector<vComplexF> vAlgebraVectorF;
-  typedef iAlgebraVector<vComplexD> vAlgebraVectorD;
-
-  typedef Lattice<vMatrix> LatticeMatrix;
-  typedef Lattice<vMatrixF> LatticeMatrixF;
-  typedef Lattice<vMatrixD> LatticeMatrixD;
-  
-  typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
-  typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
-  typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
-   
-  typedef iSUnAlgebraMatrix<vComplex>  vAlgebraMatrix;
-  typedef iSUnAlgebraMatrix<vComplexF> vAlgebraMatrixF;
-  typedef iSUnAlgebraMatrix<vComplexD> vAlgebraMatrixD;
-
-  typedef Lattice<vAlgebraMatrix>  LatticeAlgebraMatrix;
-  typedef Lattice<vAlgebraMatrixF> LatticeAlgebraMatrixF;
-  typedef Lattice<vAlgebraMatrixD> LatticeAlgebraMatrixD;
-  
-
-  typedef iSU2Matrix<Complex> SU2Matrix;
-  typedef iSU2Matrix<ComplexF> SU2MatrixF;
-  typedef iSU2Matrix<ComplexD> SU2MatrixD;
-
-  typedef iSU2Matrix<vComplex> vSU2Matrix;
-  typedef iSU2Matrix<vComplexF> vSU2MatrixF;
-  typedef iSU2Matrix<vComplexD> vSU2MatrixD;
-
-  typedef Lattice<vSU2Matrix> LatticeSU2Matrix;
-  typedef Lattice<vSU2MatrixF> LatticeSU2MatrixF;
-  typedef Lattice<vSU2MatrixD> LatticeSU2MatrixD;
-
-  // Private implementation details are specified in the following files:
-  // Grid/qcd/utils/SUn.impl
-  // Grid/qcd/utils/SUn.impl
-  // The public part of the interface follows below and refers to these
-  // private member functions.
-
-#include <Grid/qcd/utils/SUn.impl.h>
-#include <Grid/qcd/utils/Sp2n.impl.h>
-
- public:
-  template <class cplx>
-  static void generator(int lieIndex, iGroupMatrix<cplx> &ta) {
-    return generator(lieIndex, ta, group_name());
-  }
-
-  static accelerator_inline void su2SubGroupIndex(int &i1, int &i2, int su2_index) {
-    return su2SubGroupIndex(i1, i2, su2_index, group_name());
-  }
-
-  static void testGenerators(void) { testGenerators(group_name()); }
-
-  static void printGenerators(void) {
-    for (int gen = 0; gen < AlgebraDimension; gen++) {
-      Matrix ta;
-      generator(gen, ta);
-      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
-                << std::endl;
-      std::cout << GridLogMessage << ta << std::endl;
-    }
-  }
-
-  template <typename LatticeMatrixType>
-  static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out,
-                           double scale = 1.0) {
-    GridBase *grid = out.Grid();
-
-    typedef typename LatticeMatrixType::vector_type vector_type;
-
-    typedef iSinglet<vector_type> vTComplexType;
-
-    typedef Lattice<vTComplexType> LatticeComplexType;
-    typedef typename GridTypeMapper<
-        typename LatticeMatrixType::vector_object>::scalar_object MatrixType;
-
-    LatticeComplexType ca(grid);
-    LatticeMatrixType lie(grid);
-    LatticeMatrixType la(grid);
-    ComplexD ci(0.0, scale);
-    MatrixType ta;
-
-    lie = Zero();
-
-    for (int a = 0; a < AlgebraDimension; a++) {
-      random(pRNG, ca);
-
-      ca = (ca + conjugate(ca)) * 0.5;
-      ca = ca - 0.5;
-
-      generator(a, ta);
-
-      la = ci * ca * ta;
-
-      lie = lie + la;  // e^{i la ta}
-    }
-    taExp(lie, out);
-  }
-
-  static void GaussianFundamentalLieAlgebraMatrix(GridParallelRNG &pRNG,
-                                                  LatticeMatrix &out,
-                                                  Real scale = 1.0) {
-    GridBase *grid = out.Grid();
-    LatticeReal ca(grid);
-    LatticeMatrix la(grid);
-    Complex ci(0.0, scale);
-    Matrix ta;
-
-    out = Zero();
-    for (int a = 0; a < AlgebraDimension; a++) {
-      gaussian(pRNG, ca);
-      generator(a, ta);
-      la = toComplex(ca) * ta;
-      out += la;
-    }
-    out *= ci;
-  }
-
-  static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h,
-                                          LatticeMatrix &out,
-                                          Real scale = 1.0) {
-    conformable(h, out);
-    GridBase *grid = out.Grid();
-    LatticeMatrix la(grid);
-    Matrix ta;
-
-    out = Zero();
-    for (int a = 0; a < AlgebraDimension; a++) {
-      generator(a, ta);
-      la = peekColour(h, a) * timesI(ta) * scale;
-      out += la;
-    }
-  }
-
-  // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1
-  // ) inverse operation: FundamentalLieAlgebraMatrix
-  static void projectOnAlgebra(LatticeAlgebraVector &h_out,
-                               const LatticeMatrix &in, Real scale = 1.0) {
-    conformable(h_out, in);
-    h_out = Zero();
-    Matrix Ta;
-
-    for (int a = 0; a < AlgebraDimension; a++) {
-      generator(a, Ta);
-      pokeColour(h_out, -2.0 * (trace(timesI(Ta) * in)) * scale, a);
-    }
-  }
-
-   
-  template <class vtype>
-  accelerator_inline static iScalar<vtype> ProjectOnGeneralGroup(const iScalar<vtype> &r) {
-    return ProjectOnGeneralGroup(r, group_name());
-  }
-
-  template <class vtype, int N>
-  accelerator_inline static iVector<vtype,N> ProjectOnGeneralGroup(const iVector<vtype,N> &r) {
-    return ProjectOnGeneralGroup(r, group_name());
-  }
-
-  template <class vtype,int N, typename std::enable_if< GridTypeMapper<vtype>::TensorLevel == 0 >::type * =nullptr>
-  accelerator_inline static iMatrix<vtype,N> ProjectOnGeneralGroup(const iMatrix<vtype,N> &arg) {
-    return ProjectOnGeneralGroup(arg, group_name());
-  }
-
-  template <int N,class vComplex_t>                  // Projects on the general groups U(N), Sp(2N)xZ2 i.e. determinant is allowed a complex phase.
-  static void ProjectOnGeneralGroup(Lattice<iVector<iScalar<iMatrix<vComplex_t, N> >, Nd> > &U) {
-    for (int mu = 0; mu < Nd; mu++) {
-      auto Umu = PeekIndex<LorentzIndex>(U, mu);
-      Umu = ProjectOnGeneralGroup(Umu);
-    }
-  }
-       
-
-  
-  template <int N,class vComplex_t>
-  static Lattice<iScalar<iScalar<iMatrix<vComplex_t, N> > > > ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplex_t, N> > > > &Umu) {
-    return ProjectOnGeneralGroup(Umu, group_name());
-  }
-
-  template <int N,class vComplex_t>       // Projects on SU(N), Sp(2N), with unit determinant, by first projecting on general group and then enforcing unit determinant
-  static void ProjectOnSpecialGroup(Lattice<iScalar<iScalar<iMatrix<vComplex_t, N> > > > &Umu) {
-       Umu = ProjectOnGeneralGroup(Umu);
-       auto det = Determinant(Umu);
-
-       det = conjugate(det);
-
-       for (int i = 0; i < N; i++) {
-           auto element = PeekIndex<ColourIndex>(Umu, N - 1, i);
-           element = element * det;
-           PokeIndex<ColourIndex>(Umu, element, Nc - 1, i);
-       }
-   }
-
-  template <int N,class vComplex_t>    // reunitarise, resimplectify... previously ProjectSUn
-    static void ProjectOnSpecialGroup(Lattice<iVector<iScalar<iMatrix<vComplex_t, N> >, Nd> > &U) {
-      // Reunitarise
-      for (int mu = 0; mu < Nd; mu++) {
-        auto Umu = PeekIndex<LorentzIndex>(U, mu);
-        ProjectOnSpecialGroup(Umu);
-        PokeIndex<LorentzIndex>(U, Umu, mu);
-      }
-    }
-    
-  template <typename GaugeField>
-  static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iGroupMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    LatticeMatrixType tmp(out.Grid());
-    for (int mu = 0; mu < Nd; mu++) {
-      //      LieRandomize(pRNG, Umu, 1.0);
-      //      PokeIndex<LorentzIndex>(out, Umu, mu);
-      gaussian(pRNG,Umu);
-      tmp = Ta(Umu);
-      taExp(tmp,Umu);
-      ProjectOnSpecialGroup(Umu);
-      //      ProjectSUn(Umu);
-      PokeIndex<LorentzIndex>(out, Umu, mu);
-    }
-  }
-  template <typename GaugeField>
-  static void TepidConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iGroupMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    for (int mu = 0; mu < Nd; mu++) {
-      LieRandomize(pRNG, Umu, 0.01);
-      PokeIndex<LorentzIndex>(out, Umu, mu);
-    }
-  }
-    
-  template <typename GaugeField>
-  static void ColdConfiguration(GaugeField &out) {
-    typedef typename GaugeField::vector_type vector_type;
-    typedef iGroupMatrix<vector_type> vMatrixType;
-    typedef Lattice<vMatrixType> LatticeMatrixType;
-
-    LatticeMatrixType Umu(out.Grid());
-    Umu = 1.0;
-    for (int mu = 0; mu < Nd; mu++) {
-      PokeIndex<LorentzIndex>(out, Umu, mu);
-    }
-  }
-    
-  template <typename GaugeField>
-  static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
-    ColdConfiguration(out);
-  }
-
-  template <typename LatticeMatrixType>
-  static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) {
-    taProj(in, out, group_name());
-  }
-    
-  template <typename LatticeMatrixType>
-  static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) {
-    typedef typename LatticeMatrixType::scalar_type ComplexType;
-
-    LatticeMatrixType xn(x.Grid());
-    RealD nfac = 1.0;
-
-    xn = x;
-    ex = xn + ComplexType(1.0);  // 1+x
-
-    // Do a 12th order exponentiation
-    for (int i = 2; i <= 12; ++i) {
-      nfac = nfac / RealD(i);  // 1/2, 1/2.3 ...
-      xn = xn * x;             // x2, x3,x4....
-      ex = ex + xn * nfac;     // x2/2!, x3/3!....
-    }
-  }
-
-// Ta are hermitian (?)
-// Anti herm is i Ta basis
-static void LieAlgebraProject(LatticeAlgebraMatrix &out,const LatticeMatrix &in, int b)
-{
-  conformable(in, out);
-  GridBase *grid = out.Grid();
-  LatticeComplex tmp(grid);
-  Matrix ta;
-  // Using Luchang's projection convention
-  //  2 Tr{Ta Tb} A_b= 2/2 delta ab A_b = A_a
-  autoView(out_v,out,AcceleratorWrite);
-  autoView(in_v,in,AcceleratorRead);
-  int N = ncolour;
-  int NNm1 = N * (N - 1);
-  int hNNm1= NNm1/2;
-  RealD sqrt_2 = sqrt(2.0);
-  Complex ci(0.0,1.0);
-  for(int su2Index=0;su2Index<hNNm1;su2Index++){
-    int i1, i2;
-    su2SubGroupIndex(i1, i2, su2Index);
-    int ax = su2Index*2;
-    int ay = su2Index*2+1;
-    accelerator_for(ss,grid->oSites(),1,{
-	// in is traceless ANTI-hermitian whereas Grid generators are Hermitian.
-	// trace( Ta x Ci in)
-	// Bet I need to move to real part with mult by -i
-	out_v[ss]()()(ax,b) = 0.5*(real(in_v[ss]()()(i2,i1)) - real(in_v[ss]()()(i1,i2)));
-	out_v[ss]()()(ay,b) = 0.5*(imag(in_v[ss]()()(i1,i2)) + imag(in_v[ss]()()(i2,i1)));
-      });
-  }
-  for(int diagIndex=0;diagIndex<N-1;diagIndex++){
-    int k = diagIndex + 1; // diagIndex starts from 0
-    int a = NNm1+diagIndex;
-    RealD scale = 1.0/sqrt(2.0*k*(k+1));
-    accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
-	auto tmp = in_v[ss]()()(0,0);
-	for(int i=1;i<k;i++){
-	  tmp=tmp+in_v[ss]()()(i,i);
-	}
-	tmp = tmp - in_v[ss]()()(k,k)*k;
-	out_v[ss]()()(a,b) =imag(tmp) * scale;
-      });
-    }
-}
-
-  
-};
-    
-template <int ncolour>
-using SU = GaugeGroup<ncolour, GroupName::SU>;
-
-template <int ncolour>
-using Sp = GaugeGroup<ncolour, GroupName::Sp>;
-
-typedef SU<2> SU2;
-typedef SU<3> SU3;
-typedef SU<4> SU4;
-typedef SU<5> SU5;
-
-typedef SU<Nc> FundamentalMatrices;
-    
-typedef Sp<2> Sp2;
-typedef Sp<4> Sp4;
-typedef Sp<6> Sp6;
-typedef Sp<8> Sp8;
-
-template <int N,class vComplex_t>
-static void ProjectSUn(Lattice<iScalar<iScalar<iMatrix<vComplex_t, N> > > > &Umu)
-{
-    GaugeGroup<N,GroupName::SU>::ProjectOnSpecialGroup(Umu);
-}
-  
-template <int N,class vComplex_t>
-static void ProjectSUn(Lattice<iVector<iScalar<iMatrix<vComplex_t, N> >,Nd> > &U)
-{
-    GaugeGroup<N,GroupName::SU>::ProjectOnSpecialGroup(U);
-}
-    
-template <int N,class vComplex_t>
-static void ProjectSpn(Lattice<iScalar<iScalar<iMatrix<vComplex_t, N> > > > &Umu)
-{
-    GaugeGroup<N,GroupName::Sp>::ProjectOnSpecialGroup(Umu);
-}
-    
-template <int N,class vComplex_t>
-static void ProjectSpn(Lattice<iVector<iScalar<iMatrix<vComplex_t, N> >,Nd> > &U)
-{
-    GaugeGroup<N,GroupName::Sp>::ProjectOnSpecialGroup(U);
-}
-
-// Explicit specialisation for SU(3).
-static void ProjectSU3(Lattice<iScalar<iScalar<iMatrix<vComplexD, 3> > > > &Umu)
-{
-  GridBase *grid = Umu.Grid();
-  const int x = 0;
-  const int y = 1;
-  const int z = 2;
-  // Reunitarise
-  Umu = ProjectOnGroup(Umu);
-  autoView(Umu_v, Umu, CpuWrite);
-  thread_for(ss, grid->oSites(), {
-    auto cm = Umu_v[ss];
-    cm()()(2, x) = adj(cm()()(0, y) * cm()()(1, z) -
-                       cm()()(0, z) * cm()()(1, y));  // x= yz-zy
-    cm()()(2, y) = adj(cm()()(0, z) * cm()()(1, x) -
-                       cm()()(0, x) * cm()()(1, z));  // y= zx-xz
-    cm()()(2, z) = adj(cm()()(0, x) * cm()()(1, y) -
-                       cm()()(0, y) * cm()()(1, x));  // z= xy-yx
-    Umu_v[ss] = cm;
-  });
-}
-static void ProjectSU3(Lattice<iVector<iScalar<iMatrix<vComplexD, 3> >, Nd> > &U)
-{
-  GridBase *grid = U.Grid();
-  // Reunitarise
-  for (int mu = 0; mu < Nd; mu++) {
-    auto Umu = PeekIndex<LorentzIndex>(U, mu);
-    Umu = ProjectOnGroup(Umu);
-    ProjectSU3(Umu);
-    PokeIndex<LorentzIndex>(U, Umu, mu);
-  }
-}
-
-NAMESPACE_END(Grid);
-#endif

Grid/qcd/utils/GaugeGroupTwoIndex.h

@@ -1,371 +0,0 @@
-////////////////////////////////////////////////////////////////////////
-//
-// * Two index representation generators
-//
-// * Normalisation for the fundamental generators:
-//   trace ta tb = 1/2 delta_ab = T_F delta_ab
-//   T_F = 1/2  for SU(N) groups
-//
-//
-//   base for NxN two index (anti-symmetric) matrices
-//   normalized to 1 (d_ij is the kroenecker delta)
-//
-//   (e^(ij)_{kl} = 1 / sqrt(2) (d_ik d_jl +/- d_jk d_il)
-//
-//   Then the generators are written as
-//
-//   (iT_a)^(ij)(lk) = i * ( tr[e^(ij)^dag e^(lk) T^trasp_a] +
-//   tr[e^(lk)e^(ij)^dag T_a] )  //
-//
-//
-////////////////////////////////////////////////////////////////////////
-
-// Authors: David Preti, Guido Cossu
-
-#ifndef QCD_UTIL_GAUGEGROUPTWOINDEX_H
-#define QCD_UTIL_GAUGEGROUPTWOINDEX_H
-
-NAMESPACE_BEGIN(Grid);
-
-enum TwoIndexSymmetry { Symmetric = 1, AntiSymmetric = -1 };
-
-constexpr inline Real delta(int a, int b) { return (a == b) ? 1.0 : 0.0; }
-
-namespace detail {
-
-template <class cplx, int nc, TwoIndexSymmetry S>
-struct baseOffDiagonalSpHelper;
-
-template <class cplx, int nc>
-struct baseOffDiagonalSpHelper<cplx, nc, AntiSymmetric> {
-  static const int ngroup = nc / 2;
-  static void baseOffDiagonalSp(int i, int j, iScalar<iScalar<iMatrix<cplx, nc> > > &eij) {
-    eij = Zero();
-    RealD tmp;
-
-    if ((i == ngroup + j) && (1 <= j) && (j < ngroup)) {
-      for (int k = 0; k < j+1; k++) {
-        if (k < j) {
-          tmp = 1 / sqrt(j * (j + 1));
-          eij()()(k, k + ngroup) = tmp;
-          eij()()(k + ngroup, k) = -tmp;
-        }
-        if (k == j) {
-          tmp = -j / sqrt(j * (j + 1));
-          eij()()(k, k + ngroup) = tmp;
-          eij()()(k + ngroup, k) = -tmp;
-        }
-      }
-
-    }
-
-    else if (i != ngroup + j) {
-      for (int k = 0; k < nc; k++)
-        for (int l = 0; l < nc; l++) {
-          eij()()(l, k) =
-              delta(i, k) * delta(j, l) - delta(j, k) * delta(i, l);
-        }
-    }
-    RealD nrm = 1. / std::sqrt(2.0);
-    eij = eij * nrm;
-  }
-};
-
-template <class cplx, int nc>
-struct baseOffDiagonalSpHelper<cplx, nc, Symmetric> {
-  static void baseOffDiagonalSp(int i, int j, iScalar<iScalar<iMatrix<cplx, nc> > > &eij) {
-    eij = Zero();
-    for (int k = 0; k < nc; k++)
-      for (int l = 0; l < nc; l++)
-        eij()()(l, k) =
-            delta(i, k) * delta(j, l) + delta(j, k) * delta(i, l);
-
-    RealD nrm = 1. / std::sqrt(2.0);
-    eij = eij * nrm;
-  }
-};
-
-}   // closing detail namespace
-
-template <int ncolour, TwoIndexSymmetry S, class group_name>
-class GaugeGroupTwoIndex : public GaugeGroup<ncolour, group_name> {
- public:
-  // The chosen convention is that we are taking ncolour to be N in SU<N> but 2N
-  // in Sp(2N). ngroup is equal to N for SU but 2N/2 = N for Sp(2N).
-  static_assert(std::is_same<group_name, GroupName::SU>::value or
-                    std::is_same<group_name, GroupName::Sp>::value,
-                "ngroup is only implemented for SU and Sp currently.");
-  static const int ngroup =
-      std::is_same<group_name, GroupName::SU>::value ? ncolour : ncolour / 2;
-  static const int Dimension =
-      (ncolour * (ncolour + S) / 2) + (std::is_same<group_name, GroupName::Sp>::value ? (S - 1) / 2 : 0);
-  static const int DimensionAS =
-      (ncolour * (ncolour - 1) / 2) + (std::is_same<group_name, GroupName::Sp>::value ? (- 1) : 0);
-  static const int DimensionS =
-      ncolour * (ncolour + 1) / 2;
-  static const int NumGenerators =
-      GaugeGroup<ncolour, group_name>::AlgebraDimension;
-
-  template <typename vtype>
-  using iGroupTwoIndexMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
-
-  typedef iGroupTwoIndexMatrix<Complex> TIMatrix;
-  typedef iGroupTwoIndexMatrix<ComplexF> TIMatrixF;
-  typedef iGroupTwoIndexMatrix<ComplexD> TIMatrixD;
-
-  typedef iGroupTwoIndexMatrix<vComplex> vTIMatrix;
-  typedef iGroupTwoIndexMatrix<vComplexF> vTIMatrixF;
-  typedef iGroupTwoIndexMatrix<vComplexD> vTIMatrixD;
-
-  typedef Lattice<vTIMatrix> LatticeTwoIndexMatrix;
-  typedef Lattice<vTIMatrixF> LatticeTwoIndexMatrixF;
-  typedef Lattice<vTIMatrixD> LatticeTwoIndexMatrixD;
-
-  typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> >
-      LatticeTwoIndexField;
-  typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> >
-      LatticeTwoIndexFieldF;
-  typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> >
-      LatticeTwoIndexFieldD;
-
-  template <typename vtype>
-  using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
-
-  typedef iGroupMatrix<Complex> Matrix;
-  typedef iGroupMatrix<ComplexF> MatrixF;
-  typedef iGroupMatrix<ComplexD> MatrixD;
-    
-private:
-  template <class cplx>
-  static void baseDiagonal(int Index, iGroupMatrix<cplx> &eij) {
-    eij = Zero();
-    eij()()(Index - ncolour * (ncolour - 1) / 2,
-            Index - ncolour * (ncolour - 1) / 2) = 1.0;
-  }
-    
-  template <class cplx>
-  static void baseOffDiagonal(int i, int j, iGroupMatrix<cplx> &eij, GroupName::SU) {
-    eij = Zero();
-    for (int k = 0; k < ncolour; k++)
-      for (int l = 0; l < ncolour; l++)
-        eij()()(l, k) =
-            delta(i, k) * delta(j, l) + S * delta(j, k) * delta(i, l);
-
-    RealD nrm = 1. / std::sqrt(2.0);
-    eij = eij * nrm;
-  }
-    
-  template <class cplx>
-  static void baseOffDiagonal(int i, int j, iGroupMatrix<cplx> &eij, GroupName::Sp) {
-    detail::baseOffDiagonalSpHelper<cplx, ncolour, S>::baseOffDiagonalSp(i, j, eij);
-  }
-
-public:
-    
-  template <class cplx>
-  static void base(int Index, iGroupMatrix<cplx> &eij) {
-  // returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
-    assert(Index < Dimension);
-    eij = Zero();
-  // for the linearisation of the 2 indexes
-    static int a[ncolour * (ncolour - 1) / 2][2];  // store the a <-> i,j
-    static bool filled = false;
-    if (!filled) {
-      int counter = 0;
-      for (int i = 1; i < ncolour; i++) {
-      for (int j = 0; j < i; j++) {
-        if (std::is_same<group_name, GroupName::Sp>::value)
-          {
-            if (j==0 && i==ngroup+j && S==-1) {
-            //std::cout << "skipping" << std::endl; // for Sp2n this vanishes identically.
-              j = j+1;
-            }
-          }
-          a[counter][0] = i;
-          a[counter][1] = j;
-          counter++;
-          }
-      }
-      filled = true;
-    }
-    if (Index < ncolour*ncolour - DimensionS)
-    {
-      baseOffDiagonal(a[Index][0], a[Index][1], eij, group_name());
-    } else {
-      baseDiagonal(Index, eij);
-    }
-  }
-    
-  static void printBase(void) {
-    for (int gen = 0; gen < Dimension; gen++) {
-      Matrix tmp;
-      base(gen, tmp);
-      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
-                << std::endl;
-      std::cout << GridLogMessage << tmp << std::endl;
-    }
-  }
-
-  template <class cplx>
-  static void generator(int Index, iGroupTwoIndexMatrix<cplx> &i2indTa) {
-    Vector<iGroupMatrix<cplx> > ta(NumGenerators);
-    Vector<iGroupMatrix<cplx> > eij(Dimension);
-    iGroupMatrix<cplx> tmp;
-
-    for (int a = 0; a < NumGenerators; a++)
-      GaugeGroup<ncolour, group_name>::generator(a, ta[a]);
-
-    for (int a = 0; a < Dimension; a++) base(a, eij[a]);
-
-    for (int a = 0; a < Dimension; a++) {
-      tmp = transpose(eij[a]*ta[Index]) + transpose(eij[a]) * ta[Index];
-      for (int b = 0; b < Dimension; b++) {
-        Complex iTr = TensorRemove(timesI(trace(tmp * eij[b])));
-        i2indTa()()(a, b) = iTr;
-      }
-    }
-  }
-
-  static void printGenerators(void) {
-    for (int gen = 0; gen < NumGenerators; gen++) {
-      TIMatrix i2indTa;
-      generator(gen, i2indTa);
-      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
-                << std::endl;
-      std::cout << GridLogMessage << i2indTa << std::endl;
-    }
-  }
-
-  static void testGenerators(void) {
-    TIMatrix i2indTa, i2indTb;
-    std::cout << GridLogMessage << "2IndexRep - Checking if traceless"
-              << std::endl;
-    for (int a = 0; a < NumGenerators; a++) {
-      generator(a, i2indTa);
-      std::cout << GridLogMessage << a << std::endl;
-      assert(norm2(trace(i2indTa)) < 1.0e-6);
-    }
-    std::cout << GridLogMessage << std::endl;
-
-    std::cout << GridLogMessage << "2IndexRep - Checking if antihermitean"
-              << std::endl;
-    for (int a = 0; a < NumGenerators; a++) {
-      generator(a, i2indTa);
-      std::cout << GridLogMessage << a << std::endl;
-      assert(norm2(adj(i2indTa) + i2indTa) < 1.0e-6);
-    }
-
-    std::cout << GridLogMessage << std::endl;
-    std::cout << GridLogMessage
-              << "2IndexRep - Checking Tr[Ta*Tb]=delta(a,b)*(N +- 2)/2"
-              << std::endl;
-    for (int a = 0; a < NumGenerators; a++) {
-      for (int b = 0; b < NumGenerators; b++) {
-        generator(a, i2indTa);
-        generator(b, i2indTb);
-
-        // generator returns iTa, so we need a minus sign here
-        Complex Tr = -TensorRemove(trace(i2indTa * i2indTb));
-        std::cout << GridLogMessage << "a=" << a << "b=" << b << "Tr=" << Tr
-                  << std::endl;
-        if (a == b) {
-          assert(real(Tr) - ((ncolour + S * 2) * 0.5) < 1e-8);
-        } else {
-          assert(real(Tr) < 1e-8);
-        }
-        assert(imag(Tr) < 1e-8);
-      }
-    }
-    std::cout << GridLogMessage << std::endl;
-  }
-
-  static void TwoIndexLieAlgebraMatrix(
-      const typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h,
-      LatticeTwoIndexMatrix &out, Real scale = 1.0) {
-    conformable(h, out);
-    GridBase *grid = out.Grid();
-    LatticeTwoIndexMatrix la(grid);
-    TIMatrix i2indTa;
-
-    out = Zero();
-    for (int a = 0; a < NumGenerators; a++) {
-      generator(a, i2indTa);
-      la = peekColour(h, a) * i2indTa;
-      out += la;
-    }
-    out *= scale;
-  }
-
-  // Projects the algebra components
-  // of a lattice matrix ( of dimension ncol*ncol -1 )
-  static void projectOnAlgebra(
-      typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h_out,
-      const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
-    conformable(h_out, in);
-    h_out = Zero();
-    TIMatrix i2indTa;
-    Real coefficient = -2.0 / (ncolour + 2 * S) * scale;
-    // 2/(Nc +/- 2) for the normalization of the trace in the two index rep
-    for (int a = 0; a < NumGenerators; a++) {
-      generator(a, i2indTa);
-      pokeColour(h_out, real(trace(i2indTa * in)) * coefficient, a);
-    }
-  }
-
-  // a projector that keeps the generators stored to avoid the overhead of
-  // recomputing them
-  static void projector(
-      typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h_out,
-      const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
-    conformable(h_out, in);
-    // to store the generators
-    static std::vector<TIMatrix> i2indTa(NumGenerators);
-    h_out = Zero();
-    static bool precalculated = false;
-    if (!precalculated) {
-      precalculated = true;
-      for (int a = 0; a < NumGenerators; a++) generator(a, i2indTa[a]);
-    }
-
-    Real coefficient =
-        -2.0 / (ncolour + 2 * S) * scale;  // 2/(Nc +/- 2) for the normalization
-    // of the trace in the two index rep
-
-    for (int a = 0; a < NumGenerators; a++) {
-      auto tmp = real(trace(i2indTa[a] * in)) * coefficient;
-      pokeColour(h_out, tmp, a);
-    }
-  }
-};
-
-template <int ncolour, TwoIndexSymmetry S>
-using SU_TwoIndex = GaugeGroupTwoIndex<ncolour, S, GroupName::SU>;
-
-// Some useful type names
-typedef SU_TwoIndex<Nc, Symmetric> TwoIndexSymmMatrices;
-typedef SU_TwoIndex<Nc, AntiSymmetric> TwoIndexAntiSymmMatrices;
-
-typedef SU_TwoIndex<2, Symmetric> SU2TwoIndexSymm;
-typedef SU_TwoIndex<3, Symmetric> SU3TwoIndexSymm;
-typedef SU_TwoIndex<4, Symmetric> SU4TwoIndexSymm;
-typedef SU_TwoIndex<5, Symmetric> SU5TwoIndexSymm;
-
-typedef SU_TwoIndex<2, AntiSymmetric> SU2TwoIndexAntiSymm;
-typedef SU_TwoIndex<3, AntiSymmetric> SU3TwoIndexAntiSymm;
-typedef SU_TwoIndex<4, AntiSymmetric> SU4TwoIndexAntiSymm;
-typedef SU_TwoIndex<5, AntiSymmetric> SU5TwoIndexAntiSymm;
-
-template <int ncolour, TwoIndexSymmetry S>
-using Sp_TwoIndex = GaugeGroupTwoIndex<ncolour, S, GroupName::Sp>;
-
-typedef Sp_TwoIndex<Nc, Symmetric> SpTwoIndexSymmMatrices;
-typedef Sp_TwoIndex<Nc, AntiSymmetric> SpTwoIndexAntiSymmMatrices;
-
-typedef Sp_TwoIndex<2, Symmetric> Sp2TwoIndexSymm;
-typedef Sp_TwoIndex<4, Symmetric> Sp4TwoIndexSymm;
-
-typedef Sp_TwoIndex<4, AntiSymmetric> Sp4TwoIndexAntiSymm;
-
-NAMESPACE_END(Grid);
-
-#endif