Fix Zero() used on thrust::complex in WordBundle4 initialisation

Grid's Zero() sentinel is not assignable to thrust::complex<double>;
use scalarD(0) instead.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

CLAUDE.md

@@ -0,0 +1,98 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## What This Is
+
+Grid is a data-parallel C++ library for lattice QCD. It provides SIMD-vectorised lattice containers, MPI-based domain decomposition, GPU acceleration (CUDA/HIP/SYCL), and a full suite of QCD algorithms including HMC.
+
+## Build
+
+Uses GNU Autotools. The bootstrap step only needs to run once (or after `configure.ac` changes).
+
+```bash
+./bootstrap.sh                    # downloads Eigen 3.4.0, generates configure
+mkdir build && cd build
+../configure [options]
+make -j$(nproc)
+make check                        # run root-level tests
+make install
+```
+
+Key configure options:
+
+| Option | Common values |
+|--------|---------------|
+| `--enable-simd=` | `AVX2`, `AVX512`, `KNL`, `A64FX`, `NEONv8`, `GPU` |
+| `--enable-comms=` | `mpi-auto`, `mpi3-auto`, `none` |
+| `--enable-accelerator=` | `cuda`, `hip`, `sycl` |
+| `--enable-shm=` | `shmopen`, `hugetlbfs`, `nvlink` |
+| `--enable-Nc=` | `3` (default), `2`, `4`, `5` |
+| `--with-gmp=`, `--with-mpfr=`, `--with-fftw=`, `--with-lime=` | paths to libs |
+| `--enable-hdf5`, `--enable-mkl`, `--enable-lapack` | optional features |
+
+Platform recipes from `README.md`:
+- **KNL**: `--enable-simd=KNL --enable-comms=mpi3-auto --enable-mkl`
+- **Skylake/Haswell**: `--enable-simd=AVX512` or `AVX2` + `--enable-comms=mpi3-auto`
+- **AMD EPYC**: `--enable-simd=AVX2 --enable-comms=mpi3`
+- **A64FX (Fugaku)**: `--enable-simd=A64FX --enable-comms=mpi3 --enable-shm=shmget` (see `SVE_README.txt`)
+
+Required external libs: GMP, MPFR, OpenSSL, zlib.
+
+## Running Tests
+
+```bash
+# From build directory
+make check                          # root-level tests (Test_simd, Test_cshift, etc.)
+make -C tests/<subdir> tests        # build tests in a subdirectory
+./tests/core/Test_simd              # run a single test binary directly
+```
+
+Test subdirectories and their focus: `core` (SIMD, stencil, comms), `solver` (CG, GMRES, eigensolvers), `hmc` (MD integrators), `forces` (fermion forces), `lanczos`, `IO`, `smearing`, `sp2n`, `debug`.
+
+## Architecture
+
+### Layer stack (bottom to top)
+
+1. **SIMD layer** (`Grid/simd/`) — platform-specific intrinsics wrapped into `vRealF`, `vComplexD`, etc. The SIMD width and layout are compile-time constants controlled by `--enable-simd`.
+
+2. **Tensor layer** (`Grid/tensors/`) — Lorentz/colour/spin tensor algebra built on top of SIMD types. `iMatrix`, `iVector`, `iScalar` templates compose into QCD types like `ColourMatrix`, `SpinColourVector`.
+
+3. **Lattice layer** (`Grid/lattice/`) — `Lattice<T>` container: a site-local tensor replicated across a distributed Cartesian grid. All arithmetic is site-parallel and expression-template-fused.
+
+4. **Cartesian/comms layer** (`Grid/cartesian/`, `Grid/communicator/`) — `GridCartesian` holds the MPI topology and local/global geometry. `Grid/cshift/` implements nearest-neighbour halo exchange; `Grid/stencil/` is the optimised multi-hop stencil used by Dirac operators.
+
+5. **Algorithm layer** (`Grid/algorithms/`) — iterative solvers (CG, GMRES, BiCGSTAB, mixed-precision), eigensolvers (Lanczos, LAPACK), FFT, smearing.
+
+6. **QCD layer** (`Grid/qcd/`) — gauge and fermion actions, HMC integrators, observables.
+
+### QCD subsystem (`Grid/qcd/`)
+
+- `action/fermion/` — Wilson, Clover, DWF (Mobius), Staggered, twisted-mass, G-parity variants
+- `action/gauge/` — Wilson gauge, Symanzik, Iwasaki, DBW2, plaquette+rect
+- `representations/` — Fundamental, Adjoint, Two-index, Sp(2n)
+- `hmc/` — Leapfrog, OMF2/OMF4 integrators; pseudofermion refreshment; Metropolis accept/reject
+- `smearing/` — APE, Stout, HEX, gradient flow
+- `observables/` — Polyakov loop, plaquette, topological charge
+
+### GPU acceleration
+
+GPU support is injected via macros (`accelerator_for`, `accelerator_for2dNB`). The `Grid/simd/` SIMD types map to scalar on GPU device code; host code paths remain vectorised. Unified virtual memory is on by default (`--enable-unified=yes`); device-aware MPI (`--enable-accelerator-aware-mpi`) avoids device→host copies on transfers.
+
+### Memory and I/O
+
+- `Grid/allocator/` — aligned/NUMA-aware allocators; caching allocator via `--enable-alloc-cache`
+- `Grid/parallelIO/` — distributed parallel reader/writer for ILDG (via LIME), SciDAC, and native binary formats
+- `Grid/serialisation/` — text, binary, HDF5, XML/JSON serialisation of arbitrary Grid objects
+
+### HMC applications
+
+`HMC/` contains production-ready HMC driver programmes (e.g. `Mobius2p1f.cc`, `DWF_plus_DSDR_nf2plus1_Shamir_Gparity.cc`). These are built separately from the library tests.
+
+## Key Conventions
+
+- **C++17** is required throughout.
+- Template structure: most classes are templated on `<_FImpl>` (fermion impl) or `<Gimpl>` (gauge impl), which encode the representation and precision. Instantiation is controlled by `--enable-fermion-instantiations`.
+- The `RealD`/`RealF`/`ComplexD`/`ComplexF` typedefs are used everywhere; avoid raw `double`/`float`.
+- Logging uses `Grid_log`, `Grid_error` macros (from `Grid/log/`); performance-critical paths use the `GRID_TRACE` / timer macros from `Grid/perfmon/`.
+- Reductions across MPI ranks go through `GridBase::GlobalSum` / `GlobalMax`; never reduce with bare MPI calls inside library code.

Grid/lattice/Lattice_reduction.h

@@ -31,6 +31,9 @@ Author: Christoph Lehner <christoph@lhnr.de>
 #if defined(GRID_SYCL)
 #include <Grid/lattice/Lattice_reduction_sycl.h>
 #endif
+#if defined(GRID_CUDA)||defined(GRID_HIP)||defined(GRID_SYCL)
+#include <Grid/lattice/Lattice_reduction_gpu_cub.h>
+#endif
 #include <Grid/lattice/Lattice_slicesum_core.h>
 
 NAMESPACE_BEGIN(Grid);

Grid/lattice/Lattice_reduction_gpu.h

@@ -198,7 +198,7 @@ __global__ void reduceKernel(const vobj *lat, sobj *buffer, Iterator n) {
 // Possibly promote to double and sum
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   typedef decltype(lat) Iterator;
@@ -224,7 +224,7 @@ inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osi
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::vector_type  vector;
   typedef typename vobj::scalar_typeD scalarD;
@@ -244,13 +244,13 @@ inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osi
 	buf[ss] = dat[ss*words+w];
       });
       
-    ret_p[w] = sumD_gpu_small(tbuf,osites);
+    ret_p[w] = sumD_gpu_small_old(tbuf,osites);
   }
   return ret;
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_objectD sobj;
   sobj ret;
@@ -261,9 +261,9 @@ inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
   int ok = getNumBlocksAndThreads(size, sizeof(sobj), numThreads, numBlocks);
   
   if ( ok ) {
-    ret = sumD_gpu_small(lat,osites);
+    ret = sumD_gpu_small_old(lat,osites);
   } else {
-    ret = sumD_gpu_large(lat,osites);
+    ret = sumD_gpu_large_old(lat,osites);
   }
   return ret;
 }
@@ -272,20 +272,20 @@ inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu(lat,osites);
+  result = sumD_gpu_old(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_large(lat,osites);
+  result = sumD_gpu_large_old(lat,osites);
   return result;
 }
 

Grid/lattice/Lattice_reduction_gpu_cub.h

@@ -0,0 +1,317 @@
+/*************************************************************************************
+    Grid physics library, www.github.com/paboyle/Grid
+    Source file: ./Grid/lattice/Lattice_reduction_gpu_cub.h
+    Copyright (C) 2015-2024
+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
+
+#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);
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Unified lattice reduction using CUB (CUDA/HIP) and sycl::reduction (SYCL).
+//
+// CUDA/HIP: one accelerator_for pass per site to extract SIMD lanes and promote to sobjD,
+// then CUB/hipCUB DeviceReduce::Reduce over the resulting array.
+//
+// rocPRIM's DeviceReduce requires warpSize(64) threads per block, each holding one element
+// in shared memory: sizeof(T)*64 must fit in sharedMemPerBlock.  Large QCD objects such as
+// LatticePropagator (sobjD = 2304 bytes, 64*2304 = 147 KB) exceed this budget.
+//
+// For those types sumD_gpu_large groups the vobj's vector_type words in bundles of 4,
+// reducing each bundle as a WordBundle4<scalarD> (64 bytes, 64*64 = 4 KB — always safe).
+// Words that do not fill a complete bundle are zero-padded.
+//
+// SYCL: sycl::reduction handles any type size through the runtime, so one path suffices.
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#if defined(GRID_CUDA) || defined(GRID_HIP)
+
+// Bundles 4 scalar_typeD values for the radix-4 large-type reduction path.
+// sizeof = 4 * sizeof(scalarD) <= 64 bytes; 64 * 64 = 4096 bytes, safely within
+// rocPRIM's shared-memory budget on all supported devices.
+template<class scalarD>
+struct WordBundle4 {
+  scalarD w[4];
+  accelerator_inline WordBundle4 operator+(const WordBundle4 &rhs) const {
+    WordBundle4 r;
+    r.w[0] = w[0] + rhs.w[0];
+    r.w[1] = w[1] + rhs.w[1];
+    r.w[2] = w[2] + rhs.w[2];
+    r.w[3] = w[3] + rhs.w[3];
+    return r;
+  }
+};
+
+// Direct CUB reduction on the full scalar_objectD.
+// Only safe when sizeof(sobjD)*64 <= device sharedMemPerBlock.
+// Do not call directly for large composite types (e.g. LatticePropagator).
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_direct(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_objectD sobjD;
+
+  deviceVector<sobjD> per_site(osites);
+  sobjD *per_site_p = &per_site[0];
+
+  accelerator_for(ss, osites, 1, {
+    sobj tmp = Reduce(lat[ss]);
+    sobjD tmpD; tmpD = tmp;
+    per_site_p[ss] = tmpD;
+  });
+
+  sobjD zero; zeroit(zero);
+  sobjD *d_out = static_cast<sobjD *>(acceleratorAllocDevice(sizeof(sobjD)));
+  void  *d_temp = nullptr;
+  size_t temp_bytes = 0;
+
+  gpuError_t gpuErr;
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, per_site_p, d_out,
+                                        (int)osites, gpucub::Sum(), zero, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_direct: DeviceReduce size query failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  d_temp = acceleratorAllocDevice(temp_bytes);
+
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, per_site_p, d_out,
+                                        (int)osites, gpucub::Sum(), zero, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_direct: DeviceReduce failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  accelerator_barrier();
+
+  sobjD result;
+  acceleratorCopyFromDevice(d_out, &result, sizeof(sobjD));
+  acceleratorFreeDevice(d_temp);
+  acceleratorFreeDevice(d_out);
+  return result;
+}
+
+// Radix-4 word-bundle path for types too large for the direct CUB path.
+// Treats vobj as words of vector_type; groups them in bundles of 4 and reduces
+// each bundle as a WordBundle4<scalarD>.  If words % 4 != 0, the final partial
+// bundle is zero-padded so all unused slots contribute zero to the sum.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::vector_type    vector;
+  typedef typename vobj::scalar_typeD   scalarD;
+  typedef typename vobj::scalar_objectD sobjD;
+  using R4 = WordBundle4<scalarD>;
+
+  const int words = sizeof(vobj) / sizeof(vector);
+  const int nfull = words / 4;
+  const int rem   = words % 4;
+
+  sobjD ret; zeroit(ret);
+  scalarD *ret_p = (scalarD *)&ret;
+
+  iScalar<vector> *idat = (iScalar<vector> *)lat;
+  deviceVector<R4> buf(osites);
+  R4 *buf_p = &buf[0];
+
+  R4 zero4;
+  zero4.w[0] = zero4.w[1] = zero4.w[2] = zero4.w[3] = scalarD(0);
+
+  R4 *d_out = static_cast<R4 *>(acceleratorAllocDevice(sizeof(R4)));
+  void  *d_temp = nullptr;
+  size_t temp_bytes = 0;
+
+  // Probe workspace size once — type R4 and count osites are fixed across all groups.
+  gpuError_t gpuErr;
+  gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                        (int)osites, gpucub::Sum(), zero4, computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "sumD_gpu_large: DeviceReduce size query failed: "
+              << gpuErr << std::endl;
+    exit(EXIT_FAILURE);
+  }
+  d_temp = acceleratorAllocDevice(temp_bytes);
+
+  // Full groups of 4 words.
+  for (int g = 0; g < nfull; g++) {
+    int base = 4 * g;
+    accelerator_for(ss, osites, 1, {
+      R4 r4;
+      r4.w[0] = TensorRemove(Reduce(idat[ss * words + base    ]));
+      r4.w[1] = TensorRemove(Reduce(idat[ss * words + base + 1]));
+      r4.w[2] = TensorRemove(Reduce(idat[ss * words + base + 2]));
+      r4.w[3] = TensorRemove(Reduce(idat[ss * words + base + 3]));
+      buf_p[ss] = r4;
+    });
+    gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                          (int)osites, gpucub::Sum(), zero4, computeStream);
+    if (gpuErr != gpuSuccess) {
+      std::cout << GridLogError << "sumD_gpu_large: DeviceReduce failed (group "
+                << g << "): " << gpuErr << std::endl;
+      exit(EXIT_FAILURE);
+    }
+    accelerator_barrier();
+    R4 group_result;
+    acceleratorCopyFromDevice(d_out, &group_result, sizeof(R4));
+    ret_p[base    ] = group_result.w[0];
+    ret_p[base + 1] = group_result.w[1];
+    ret_p[base + 2] = group_result.w[2];
+    ret_p[base + 3] = group_result.w[3];
+  }
+
+  // Partial last group: zero-pad unused slots so they contribute nothing to the sum.
+  if (rem > 0) {
+    int base = 4 * nfull;
+    accelerator_for(ss, osites, 1, {
+      R4 r4;
+      r4.w[0] = r4.w[1] = r4.w[2] = r4.w[3] = scalarD(0);
+      for (int k = 0; k < rem; k++)
+        r4.w[k] = TensorRemove(Reduce(idat[ss * words + base + k]));
+      buf_p[ss] = r4;
+    });
+    gpuErr = gpucub::DeviceReduce::Reduce(d_temp, temp_bytes, buf_p, d_out,
+                                          (int)osites, gpucub::Sum(), zero4, computeStream);
+    if (gpuErr != gpuSuccess) {
+      std::cout << GridLogError << "sumD_gpu_large: DeviceReduce failed (partial group): "
+                << gpuErr << std::endl;
+      exit(EXIT_FAILURE);
+    }
+    accelerator_barrier();
+    R4 partial_result;
+    acceleratorCopyFromDevice(d_out, &partial_result, sizeof(R4));
+    for (int k = 0; k < rem; k++)
+      ret_p[4 * nfull + k] = partial_result.w[k];
+  }
+
+  acceleratorFreeDevice(d_temp);
+  acceleratorFreeDevice(d_out);
+  return ret;
+}
+
+// Dispatch: direct CUB path for types that fit in the shared-memory budget,
+// radix-4 word-bundle path for larger types.
+// Threshold 512 bytes: 64 * 512 = 32768 bytes, within rocPRIM's
+// ROCPRIM_SHARED_MEMORY_MAX on all supported devices.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_objectD sobjD;
+  if constexpr (sizeof(sobjD) > 512) {
+    return sumD_gpu_large(lat, osites);
+  } else {
+    return sumD_gpu_direct(lat, osites);
+  }
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu(lat, osites);
+  return result;
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu_large(lat, osites);
+  return result;
+}
+
+#endif // GRID_CUDA || GRID_HIP
+
+#if defined(GRID_SYCL)
+
+// Accumulates in sobjD throughout, fixing the precision bug in the original
+// Lattice_reduction_sycl.h which accumulated in sobj then converted at the end.
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_objectD sobjD;
+
+  sobjD identity; zeroit(identity);
+  sobjD ret;      zeroit(ret);
+  {
+    sycl::buffer<sobjD, 1> abuff(&ret, {1});
+    theGridAccelerator->submit([&](sycl::handler &cgh) {
+      auto Reduction = sycl::reduction(abuff, cgh, identity, std::plus<>());
+      cgh.parallel_for(sycl::range<1>{(size_t)osites},
+                       Reduction,
+                       [=](sycl::id<1> item, auto &sum) {
+                         sobj s = Reduce(lat[item[0]]);
+                         sobjD sd; sd = s;
+                         sum += sd;
+                       });
+    });
+  }
+  return ret;
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+{
+  return sumD_gpu(lat, osites);
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites)
+{
+  typedef typename vobj::scalar_object sobj;
+  sobj result;
+  result = sumD_gpu(lat, osites);
+  return result;
+}
+
+template<class vobj>
+inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+{
+  return sum_gpu(lat, osites);
+}
+
+#endif // GRID_SYCL
+
+NAMESPACE_END(Grid);

Grid/lattice/Lattice_reduction_sycl.h

@@ -6,7 +6,7 @@ NAMESPACE_BEGIN(Grid);
 
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_tensor(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_objectD sumD_gpu_tensor_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   typedef typename vobj::scalar_objectD sobjD;
@@ -31,40 +31,40 @@ inline typename vobj::scalar_objectD sumD_gpu_tensor(const vobj *lat, Integer os
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_large_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_tensor(lat,osites);
+  return sumD_gpu_tensor_old(lat,osites);
 }
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu_small(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_small_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large(lat,osites);
+  return sumD_gpu_large_old(lat,osites);
 }
 
 template <class vobj>
-inline typename vobj::scalar_objectD sumD_gpu(const vobj *lat, Integer osites)
+inline typename vobj::scalar_objectD sumD_gpu_old(const vobj *lat, Integer osites)
 {
-  return sumD_gpu_large(lat,osites);
+  return sumD_gpu_large_old(lat,osites);
 }
 
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Return as same precision as input performing reduction in double precision though
 /////////////////////////////////////////////////////////////////////////////////////////////////////////
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu(const vobj *lat, Integer osites) 
+inline typename vobj::scalar_object sum_gpu_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu(lat,osites);
+  result = sumD_gpu_old(lat,osites);
   return result;
 }
 
 template <class vobj>
-inline typename vobj::scalar_object sum_gpu_large(const vobj *lat, Integer osites)
+inline typename vobj::scalar_object sum_gpu_large_old(const vobj *lat, Integer osites)
 {
   typedef typename vobj::scalar_object sobj;
   sobj result;
-  result = sumD_gpu_large(lat,osites);
+  result = sumD_gpu_large_old(lat,osites);
   return result;
 }
 

tests/debug/Test_reduction.cc

@@ -0,0 +1,160 @@
+/*************************************************************************************
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./tests/debug/Test_reduction.cc
+
+    Copyright (C) 2024
+
+Author: Peter Boyle <pboyle@bnl.gov>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+
+static int passed = 0;
+static int failed = 0;
+
+static void check(bool ok, const std::string &msg)
+{
+  if (ok) {
+    std::cout << GridLogMessage << "PASS  " << msg << std::endl;
+    passed++;
+  } else {
+    std::cout << GridLogMessage << "FAIL  " << msg << std::endl;
+    failed++;
+  }
+}
+
+// Frobenius-like norm squared of any Grid scalar tensor object.
+// For iScalar chains:  real(conj(a)*a)
+// For iMatrix<T,N>:    sum_{i,j} real(conj(a_ij)*a_ij)  (all elements)
+// Both old and new paths promote to double before accumulating, so the result
+// is always exact enough for the tolerances used here.
+template<class T>
+RealD scalarNorm2(const T &a)
+{
+  return (RealD)real(TensorRemove(innerProduct(a, a)));
+}
+
+template<class Field>
+void testReduction(GridCartesian *grid, GridParallelRNG &rng,
+                   const std::string &name, int Ncomp)
+{
+  typedef typename Field::vector_object vobj;
+  typedef typename vobj::scalar_object  sobj;
+  typedef typename vobj::scalar_type    scalar_type;
+
+  const Integer V      = grid->_gsites;
+  const Integer osites = grid->oSites();
+
+  // Detect single vs double precision by comparing fundamental scalar sizes.
+  const bool isFloat = (sizeof(scalar_type) < sizeof(ComplexD));
+
+  std::cout << GridLogMessage << "=== " << name << " ===" << std::endl;
+
+  Field field(grid);
+
+  //--------------------------------------------------------------------
+  // a) Gaussian random field: sum_gpu (new CUB path) vs sum_gpu_old
+  //    (preserved hand-rolled shared-memory path).  Both promote lanes
+  //    to double internally, so results should agree to near-roundoff.
+  //--------------------------------------------------------------------
+#if defined(GRID_CUDA) || defined(GRID_HIP) || defined(GRID_SYCL)
+  {
+    gaussian(rng, field);
+
+    autoView(v, field, AcceleratorRead);
+    sobj new_result = sum_gpu    (&v[0], osites);
+    sobj old_result = sum_gpu_old(&v[0], osites);
+
+    sobj  diff  = new_result - old_result;
+    RealD diffn = scalarNorm2(diff);
+    RealD refn  = scalarNorm2(old_result);
+    RealD reldiff = (refn > 0.0) ? std::sqrt(diffn / refn) : std::sqrt(diffn);
+
+    // Float fields: both paths cast from double to float, expect O(eps_float).
+    // Double fields: ordering differences at most O(V * eps_double).
+    RealD tol = isFloat ? 1e-6 : 1e-10;
+
+    std::cout << GridLogMessage
+              << name << " random reldiff = " << reldiff << std::endl;
+    check(reldiff < tol, name + " random: sum_gpu agrees with sum_gpu_old");
+  }
+#endif
+
+  //--------------------------------------------------------------------
+  // b) Constant field via field = 1.0.
+  //
+  //    Grid's iMatrix::operator=(scalar) sets only the diagonal, so:
+  //      LatticeComplex       -> scalar 1.0        (Ncomp = 1  nonzero per site)
+  //      LatticeColourMatrix  -> Nc x Nc identity  (Ncomp = Nc nonzero per site)
+  //      LatticePropagator    -> (Ns*Nc)^2 identity (Ncomp = Ns*Nc nonzero per site)
+  //
+  //    After GlobalSum: sum_result has Ncomp diagonal entries each equal to V,
+  //    all off-diagonal entries zero.  Grid's recursive innerProduct computes
+  //    the Frobenius inner product (sum of |element|^2 over all indices), giving
+  //
+  //        innerProduct(sum_result, sum_result) = Ncomp * V^2
+  //--------------------------------------------------------------------
+  {
+    field = 1.0;
+    sobj sum_result = sum(field);   // uses new GPU path + GlobalSum
+
+    RealD got      = scalarNorm2(sum_result);
+    RealD expected = (RealD)Ncomp * (RealD)V * (RealD)V;
+    RealD reldiff  = std::abs(got - expected) / expected;
+
+    std::cout << GridLogMessage
+              << name << " const: got " << got
+              << "  expected " << expected
+              << "  reldiff "  << reldiff << std::endl;
+    check(reldiff < 1e-8, name + " const: innerProduct(sum,sum) = Ncomp*V^2");
+  }
+}
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  Coordinate latt = GridDefaultLatt();
+  Coordinate simd = GridDefaultSimd(Nd, vComplexD::Nsimd());
+  Coordinate mpi  = GridDefaultMpi();
+
+  GridCartesian *grid = SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi);
+  GridParallelRNG rng(grid);
+  rng.SeedFixedIntegers({1, 2, 3, 4});
+
+  std::cout << GridLogMessage << "Lattice : " << latt << std::endl;
+  std::cout << GridLogMessage << "Volume  : " << grid->_gsites << std::endl;
+
+  testReduction<LatticeComplexF>      (grid, rng, "LatticeComplexF",      1      );
+  testReduction<LatticeComplexD>      (grid, rng, "LatticeComplexD",      1      );
+  testReduction<LatticeColourMatrixF> (grid, rng, "LatticeColourMatrixF", Nc     );
+  testReduction<LatticeColourMatrixD> (grid, rng, "LatticeColourMatrixD", Nc     );
+  testReduction<LatticePropagatorF>   (grid, rng, "LatticePropagatorF",   Ns*Nc  );
+  testReduction<LatticePropagatorD>   (grid, rng, "LatticePropagatorD",   Ns*Nc  );
+
+  std::cout << GridLogMessage << "==============================" << std::endl;
+  std::cout << GridLogMessage << passed << " PASSED   " << failed << " FAILED" << std::endl;
+
+  Grid_finalize();
+  return (failed > 0) ? EXIT_FAILURE : EXIT_SUCCESS;
+}