refactor slicesum: slicesum uses GPU version by default now

Grid/lattice/Lattice_reduction.h

@@ -27,12 +27,11 @@ Author: Christoph Lehner <christoph@lhnr.de>
 
 #if defined(GRID_CUDA)||defined(GRID_HIP)
 #include <Grid/lattice/Lattice_reduction_gpu.h>
-#include <Grid/lattice/Lattice_slicesum_gpu.h>
 #endif
 #if defined(GRID_SYCL)
 #include <Grid/lattice/Lattice_reduction_sycl.h>
-#include <Grid/lattice/Lattice_slicesum_sycl.h>
 #endif
+#include <Grid/lattice/Lattice_slicesum_core.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -450,19 +449,10 @@ 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
+  
-  // Parallel over orthog direction
+  //Reduce Data down to lvSum
-  autoView( Data_v, Data, CpuRead);
+  sliceSumReduction(Data,lvSum,rd, e1,e2,stride,ostride,Nsimd);
-  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);
@@ -506,19 +496,6 @@ sliceSum(const Lattice<vobj> &Data,int orthogdim)
   return result;
 }
 
-template<class vobj> inline
-std::vector<typename vobj::scalar_object> 
-sliceSumGpu(const Lattice<vobj> &Data,int orthogdim)
-{
-  std::vector<typename vobj::scalar_object> result;
-  #if defined(GRID_CUDA) || defined(GRID_HIP)
-  sliceSumGpu(Data,result,orthogdim);
-  #elif defined(GRID_SYCL)
-  sliceSum_sycl(Data,result,orthogdim);
-  #endif
-  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_slicesum_core.h

@@ -0,0 +1,204 @@
+#pragma once
+#if defined(GRID_CUDA)
+
+#include <cub/cub.cuh>
+#define gpucub cub
+#define gpuMalloc cudaMalloc
+#define gpuMemcpyAsync cudaMemcpyAsync
+#define gpuMemcpyDeviceToHost cudaMemcpyDeviceToHost
+#define gpuMemcpyHostToDevice cudaMemcpyHostToDevice
+#define gpuError_t cudaError_t
+#define gpuSuccess cudaSuccess
+
+#elif defined(GRID_HIP)
+
+#include <hipcub/hipcub.hpp>
+#define gpucub hipcub
+#define gpuMalloc hipMalloc
+#define gpuMemcpyAsync hipMemcpyAsync
+#define gpuMemcpyDeviceToHost hipMemcpyDeviceToHost
+#define gpuMemcpyHostToDevice hipMemcpyHostToDevice
+#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(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)
+{
+  typedef typename vobj::scalar_object sobj;
+
+  size_t subvol_size = e1*e2;
+
+  commVector<vobj> reduction_buffer(rd*subvol_size);
+  auto rb_p = &reduction_buffer[0];
+
+  vobj vobj_zero; //Need to provide initial value for reduction operation
+  zeroit(vobj_zero);
+  
+
+  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
+  gpuError_t gpuErr = gpuMalloc(&d_out,rd*sizeof(vobj));
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (d_out)! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  gpuErr = gpuMalloc(&d_offsets,sizeof(int)*(rd+1));
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (d_offsets)! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  //copy offsets to device
+  gpuErr = gpuMemcpyAsync(d_offsets,&offsets[0],sizeof(int)*(rd+1),gpuMemcpyHostToDevice,computeStream);
+  if (gpuErr != gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMemcpy (d_offsets)! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  //determine temp_storage_array size
+  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(temp_storage_array, temp_storage_bytes, rb_p,d_out, rd, d_offsets, d_offsets+1, ::gpucub::Sum(), vobj_zero, 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  
+  gpuErr = gpuMalloc(&temp_storage_array,temp_storage_bytes);
+  if (gpuErr!=gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (temp_storage_array)! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+  
+  autoView( Data_v, Data, AcceleratorRead);
+  //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_v[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(), vobj_zero, computeStream);
+  if (gpuErr!=gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+  
+  gpuErr = gpuMemcpyAsync(&lvSum[0],d_out,rd*sizeof(vobj),gpuMemcpyDeviceToHost,computeStream);
+  if (gpuErr!=gpuSuccess) {
+    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMemcpy (d_out)! Error: " << gpuErr <<std::endl;
+    exit(EXIT_FAILURE);
+  }
+
+  //sync after copy
+  accelerator_barrier();
+ 
+
+}
+#endif
+
+
+#if defined(GRID_SYCL)
+template<class vobj> inline void sliceSumReduction_sycl(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)
+{
+  typedef typename vobj::scalar_object sobj;
+  size_t subvol_size = e1*e2;
+
+  vobj *mysum = (vobj *) malloc_shared(sizeof(vobj),*theGridAccelerator);
+  vobj vobj_zero;
+  zeroit(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_v[ss]));
+
+  });
+
+  for (int r = 0; r < rd; r++) {
+      mysum[0] = vobj_zero; //dirty hack: cannot pass vobj_zero as identity to sycl::reduction as its not device_copyable
+      theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
+          auto Reduction = cl::sycl::reduction(mysum,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();
+      lvSum[r] = mysum[0];
+  }
+
+}
+#endif
+
+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)
+  
+  sliceSumReduction_cub(Data, lvSum, rd, e1, e2, stride, ostride, Nsimd);
+  
+  #elif defined(GRID_SYCL)
+  
+  sliceSumReduction_sycl(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_slicesum_gpu.h

@@ -1,180 +0,0 @@
-#pragma once
-#if defined(GRID_CUDA)
-
-#include <cub/cub.cuh>
-#define gpucub cub
-#define gpuMalloc cudaMalloc
-#define gpuMemcpyAsync cudaMemcpyAsync
-#define gpuMemcpyDeviceToHost cudaMemcpyDeviceToHost
-#define gpuMemcpyHostToDevice cudaMemcpyHostToDevice
-#define gpuError_t cudaError_t
-#define gpuSuccess cudaSuccess
-
-#elif defined(GRID_HIP)
-
-#include <hipcub/hipcub.hpp>
-#define gpucub hipcub
-#define gpuMalloc hipMalloc
-#define gpuMemcpyAsync hipMemcpyAsync
-#define gpuMemcpyDeviceToHost hipMemcpyDeviceToHost
-#define gpuMemcpyHostToDevice hipMemcpyHostToDevice
-#define gpuError_t hipError_t
-#define gpuSuccess hipSuccess
-
-#endif
-
-
-NAMESPACE_BEGIN(Grid);
-
-template<class vobj> inline void sliceSumGpu(const Lattice<vobj> &Data,std::vector<typename vobj::scalar_object> &result,int orthogdim)
-{
- 
-  typedef typename vobj::scalar_object sobj;
-  typedef typename vobj::scalar_object::scalar_type scalar_type;
-  GridBase  *grid = Data.Grid();
-  assert(grid!=NULL);
-
-  const int    Nd = grid->_ndimension;
-  const int Nsimd = grid->Nsimd();
-
-  assert(orthogdim >= 0);
-  assert(orthogdim < Nd);
-
-  int fd=grid->_fdimensions[orthogdim];
-  int ld=grid->_ldimensions[orthogdim];
-  int rd=grid->_rdimensions[orthogdim];
-
-  int e1=    grid->_slice_nblock[orthogdim];
-  int e2=    grid->_slice_block [orthogdim];
-  int stride=grid->_slice_stride[orthogdim];
-  int ostride=grid->_ostride[orthogdim];
-  size_t subvol_size = e1*e2;
-
-  Vector<vobj> lvSum(rd); 
-  Vector<sobj> lsSum(ld,Zero());                    
-  commVector<vobj> reduction_buffer(rd*e1*e2);
-  ExtractBuffer<sobj> extracted(Nsimd);                  
-  
-  result.resize(fd);
-
-  for(int r=0;r<rd;r++){
-    lvSum[r]=Zero();
-  }
-
-  vobj vobj_zero; //Need to provide initial value for reduction operation
-  zeroit(vobj_zero);
-  
-  autoView( Data_v, Data, AcceleratorRead);
-
-  auto rb_p = &reduction_buffer[0];
-  void *helperArray = NULL;
-  vobj *d_out;
-  size_t temp_storage_bytes = 0;
-  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
-  gpuError_t gpuErr = gpuMalloc(&d_out,rd*sizeof(vobj));
-  if (gpuErr != gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (d_out)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  gpuErr = gpuMalloc(&d_offsets,sizeof(int)*(rd+1));
-  if (gpuErr != gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (d_offsets)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  //copy offsets to device
-  gpuErr = gpuMemcpyAsync(d_offsets,&offsets[0],sizeof(int)*(rd+1),gpuMemcpyHostToDevice,computeStream);
-  if (gpuErr != gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMemcpy (d_offsets)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  //determine helperArray size
-  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(helperArray, temp_storage_bytes, rb_p,d_out, rd, d_offsets, d_offsets+1, ::gpucub::Sum(), vobj_zero, 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 helperArray  
-  gpuErr = gpuMalloc(&helperArray,temp_storage_bytes);
-  if (gpuErr!=gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMalloc (helperArray)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  //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, grid->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_v[ss]));
-
-  });
-  
-  //issue segmented reductions in computeStream
-  gpuErr = gpucub::DeviceSegmentedReduce::Reduce(helperArray, temp_storage_bytes, rb_p, d_out, rd, d_offsets, d_offsets+1,::gpucub::Sum(), vobj_zero, computeStream);
-  if (gpuErr!=gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpucub::DeviceSegmentedReduce::Reduce! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-  
-  gpuErr = gpuMemcpyAsync(&lvSum[0],d_out,rd*sizeof(vobj),gpuMemcpyDeviceToHost,computeStream);
-  if (gpuErr!=gpuSuccess) {
-    std::cout << GridLogError << "Lattice_slicesum_gpu.h: Encountered error during gpuMemcpy (d_out)! Error: " << gpuErr <<std::endl;
-    exit(EXIT_FAILURE);
-  }
-
-  //sync after copy
-  accelerator_barrier();
- 
-  // Sum across simd lanes in the plane, breaking out orthog dir.
-  Coordinate icoor(Nd);
-
-  for(int rt=0;rt<rd;rt++){
-
-    extract(lvSum[rt],extracted);
-
-    for(int idx=0;idx<Nsimd;idx++){
-
-      grid->iCoorFromIindex(icoor,idx);
-
-      int ldx =rt+icoor[orthogdim]*rd;
-      
-      lsSum[ldx]=lsSum[ldx]+extracted[idx];
-
-    }
-  }
-
-  // sum over nodes.
-  for(int t=0;t<fd;t++){
-    int pt = t/ld; // processor plane
-    int lt = t%ld;
-    if ( pt == grid->_processor_coor[orthogdim] ) {
-      result[t]=lsSum[lt];
-    } else {
-      result[t]=Zero();
-    }
-
-  }
-  scalar_type * ptr = (scalar_type *) &result[0];
-  int words = fd*sizeof(sobj)/sizeof(scalar_type);
-  grid->GlobalSumVector(ptr, words);
-}
-
-NAMESPACE_END(Grid);

Grid/lattice/Lattice_slicesum_sycl.h

@@ -1,110 +0,0 @@
-#pragma once
-
-NAMESPACE_BEGIN(Grid);
-
-template <class vobj>
-inline void sliceSum_sycl(const Lattice<vobj> &Data, std::vector<typename vobj::scalar_object> &result, int orthogdim)
-{
-    typedef typename vobj::scalar_object sobj;
-    typedef typename vobj::scalar_object::scalar_type scalar_type;
-
-    GridBase *grid = Data.Grid();
-    assert(grid!=NULL);
-
-    const int Nd = grid->_ndimension;
-    const size_t Nsimd = grid->Nsimd();
-
-    assert(orthogdim >= 0);
-    assert(orthogdim < Nd);
-
-    int fd=grid->_fdimensions[orthogdim];
-    int ld=grid->_ldimensions[orthogdim];
-    int rd=grid->_rdimensions[orthogdim];
-
-    int e1=    grid->_slice_nblock[orthogdim];
-    int e2=    grid->_slice_block [orthogdim];
-    int stride=grid->_slice_stride[orthogdim];
-    int ostride=grid->_ostride[orthogdim];
-    size_t subvol_size = e1*e2;
-
-    vobj *mysum = (vobj *) malloc_shared(sizeof(vobj),*theGridAccelerator);
-    
-    result.resize(fd);
-
-    Vector<vobj> lvSum(rd); 
-    Vector<sobj> lsSum(ld,Zero());                    
-    commVector<vobj> reduction_buffer(rd*subvol_size);
-    ExtractBuffer<sobj> extracted(Nsimd);      
-    vobj vobj_zero;
-    zeroit(vobj_zero);
-
-    for(int r=0;r<rd;r++){
-        lvSum[r]=Zero();
-    }
-
-    auto rb_p = &reduction_buffer[0];
-
-    autoView(Data_v, Data, AcceleratorRead);
-
-    //prepare reduction buffer 
-    accelerator_for2d( 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_v[ss]));
-
-    });
-  
-    for (int r = 0; r < rd; r++) {
-        mysum[0] = vobj_zero; //dirty hack: cannot pass vobj_zero as identity to sycl::reduction as its not device_copyable
-        theGridAccelerator->submit([&](cl::sycl::handler &cgh) {
-            auto Reduction = cl::sycl::reduction(mysum,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();
-        lvSum[r] = mysum[0];
-    }
-
-    Coordinate icoor(Nd);
-
-    for(int rt=0;rt<rd;rt++){
-
-        extract(lvSum[rt],extracted);
-
-        for(int idx=0;idx<Nsimd;idx++){
-
-        grid->iCoorFromIindex(icoor,idx);
-
-        int ldx =rt+icoor[orthogdim]*rd;
-        
-        lsSum[ldx]=lsSum[ldx]+extracted[idx];
-
-        }
-    }
-
-    // sum over nodes.
-    for(int t=0;t<fd;t++){
-        int pt = t/ld; // processor plane
-        int lt = t%ld;
-        if ( pt == grid->_processor_coor[orthogdim] ) {
-        result[t]=lsSum[lt];
-        } else {
-        result[t]=Zero();
-        }
-
-    }
-    scalar_type * ptr = (scalar_type *) &result[0];
-    int words = fd*sizeof(sobj)/sizeof(scalar_type);
-    grid->GlobalSumVector(ptr, words);
-
-}
-
-NAMESPACE_END(Grid);

tests/core/Test_sliceSum.cc

@@ -1,5 +1,79 @@
 #include <Grid/Grid.h>
 
+template<class vobj> inline void sliceSumCPU(const Grid::Lattice<vobj> &Data,std::vector<typename vobj::scalar_object> &result,int orthogdim)
+{
+  using namespace Grid;
+  ///////////////////////////////////////////////////////
+  // FIXME precision promoted summation
+  // may be important for correlation functions
+  // But easily avoided by using double precision fields
+  ///////////////////////////////////////////////////////
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_object::scalar_type scalar_type;
+  GridBase  *grid = Data.Grid();
+  assert(grid!=NULL);
+
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  assert(orthogdim >= 0);
+  assert(orthogdim < Nd);
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+
+  Vector<vobj> lvSum(rd); // will locally sum vectors first
+  Vector<sobj> lsSum(ld,Zero());                    // sum across these down to scalars
+  ExtractBuffer<sobj> extracted(Nsimd);                  // splitting the SIMD
+
+  result.resize(fd); // And then global sum to return the same vector to every node 
+  for(int r=0;r<rd;r++){
+    lvSum[r]=Zero();
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+  int ostride=grid->_ostride[orthogdim];
+  
+  //Reduce Data down to lvSum
+  sliceSumReduction_cpu(Data,lvSum,rd, e1,e2,stride,ostride,Nsimd);
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  Coordinate icoor(Nd);
+
+  for(int rt=0;rt<rd;rt++){
+
+    extract(lvSum[rt],extracted);
+
+    for(int idx=0;idx<Nsimd;idx++){
+
+      grid->iCoorFromIindex(icoor,idx);
+
+      int ldx =rt+icoor[orthogdim]*rd;
+
+      lsSum[ldx]=lsSum[ldx]+extracted[idx];
+
+    }
+  }
+  
+  // sum over nodes.
+  for(int t=0;t<fd;t++){
+    int pt = t/ld; // processor plane
+    int lt = t%ld;
+    if ( pt == grid->_processor_coor[orthogdim] ) {
+      result[t]=lsSum[lt];
+    } else {
+      result[t]=Zero();
+    }
+
+  }
+  scalar_type * ptr = (scalar_type *) &result[0];
+  int words = fd*sizeof(sobj)/sizeof(scalar_type);
+  grid->GlobalSumVector(ptr, words);
+}
+
 
 int main (int argc, char ** argv) {
     
@@ -26,7 +100,7 @@ int main (int argc, char ** argv) {
 
     //warmup
     for (int sweeps = 0; sweeps < 5; sweeps++) {
-      reduction_result = sliceSumGpu(test_data,0);
+      reduction_result = sliceSum(test_data,0);
     }
 
     int trace_id = traceStart("sliceSum benchmark");
@@ -35,23 +109,23 @@ int main (int argc, char ** argv) {
       RealD t=-usecond();
 
       tracePush("sliceSum");
-      sliceSum(test_data,reduction_reference,i);
+      sliceSumCPU(test_data,reduction_reference,i);
       tracePop("sliceSum");
 
       t+=usecond();
-
+      std::cout << GridLogMessage << "Orthog. dir. = " << i << std::endl;
-      std::cout << GridLogMessage << " sliceSum took "<<t<<" usecs"<<std::endl;
+      std::cout << GridLogMessage << "CPU sliceSum took "<<t<<" usecs"<<std::endl;
       
       
       RealD tgpu=-usecond();
 
       tracePush("sliceSumGpu");
-      reduction_result = sliceSumGpu(test_data,i);
+      reduction_result = sliceSum(test_data,i);
       tracePop("sliceSumGpu");
 
       tgpu+=usecond();
 
-      std::cout << GridLogMessage <<" sliceSumGpu took "<<tgpu<<" usecs"<<std::endl;
+      std::cout << GridLogMessage <<"GPU sliceSum took "<<tgpu<<" usecs"<<std::endl<<std::endl;;
 
 
       for(int t=0;t<reduction_reference.size();t++) {