Add an OpenCL dot kernel

We have to name the kernel stream_dot (for example) because the
"dot" kernel already exists.

OCLStream.cpp

@@ -50,6 +50,29 @@ std::string kernels{R"CLC(
     a[i] = b[i] + scalar * c[i];
   }
 
+  kernel void stream_dot(
+    global const TYPE * restrict a,
+    global const TYPE * restrict b,
+    global TYPE * restrict sum,
+    local TYPE * restrict wg_sum)
+  {
+    const size_t i = get_global_id(0);
+    const size_t local_i = get_local_id(0);
+    wg_sum[local_i] = a[i] * b[i];
+
+    for (int offset = get_local_size(0) / 2; offset > 0; offset /= 2)
+    {
+      barrier(CLK_LOCAL_MEM_FENCE);
+      if (local_i < offset)
+      {
+        wg_sum[local_i] += wg_sum[local_i+offset];
+      }
+    }
+
+    if (local_i == 0)
+      sum[get_group_id(0)] = wg_sum[local_i];
+  }
+
 )CLC"};
 
 
@@ -99,6 +122,7 @@ OCLStream<T>::OCLStream(const unsigned int ARRAY_SIZE, const int device_index)
   mul_kernel = new cl::KernelFunctor<cl::Buffer, cl::Buffer>(program, "mul");
   add_kernel = new cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer>(program, "add");
   triad_kernel = new cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer>(program, "triad");
+  dot_kernel = new cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer, cl::LocalSpaceArg>(program, "stream_dot");
 
   array_size = ARRAY_SIZE;
 
@@ -114,6 +138,7 @@ OCLStream<T>::OCLStream(const unsigned int ARRAY_SIZE, const int device_index)
   d_a = cl::Buffer(context, CL_MEM_READ_WRITE, sizeof(T) * ARRAY_SIZE);
   d_b = cl::Buffer(context, CL_MEM_READ_WRITE, sizeof(T) * ARRAY_SIZE);
   d_c = cl::Buffer(context, CL_MEM_READ_WRITE, sizeof(T) * ARRAY_SIZE);
+  d_sum = cl::Buffer(context, CL_MEM_WRITE_ONLY, sizeof(T) * WGSIZE);
 
 }
 
@@ -166,6 +191,22 @@ void OCLStream<T>::triad()
   queue.finish();
 }
 
+template <class T>
+T OCLStream<T>::dot()
+{
+  (*dot_kernel)(
+    cl::EnqueueArgs(queue, cl::NDRange(array_size), cl::NDRange(WGSIZE)),
+    d_a, d_b, d_sum, cl::Local(sizeof(T) * WGSIZE)
+  );
+  cl::copy(queue, d_sum, sums.begin(), sums.end());
+
+  T sum = 0.0;
+  for (T val : sums)
+    sum += val;
+
+  return sum;
+}
+
 template <class T>
 void OCLStream<T>::write_arrays(const std::vector<T>& a, const std::vector<T>& b, const std::vector<T>& c)
 {

OCLStream.h

@@ -20,6 +20,9 @@
 
 #define IMPLEMENTATION_STRING "OpenCL"
 
+// Local work-group size for dot kernel
+#define WGSIZE 1024
+
 template <class T>
 class OCLStream : public Stream<T>
 {
@@ -27,10 +30,14 @@ class OCLStream : public Stream<T>
     // Size of arrays
     unsigned int array_size;
     
+    // Host array for partial sums for dot kernel
+    std::vector<T> sums;
+
     // Device side pointers to arrays
     cl::Buffer d_a;
     cl::Buffer d_b;
     cl::Buffer d_c;
+    cl::Buffer d_sum;
 
     // OpenCL objects
     cl::Device device;
@@ -41,6 +48,7 @@ class OCLStream : public Stream<T>
     cl::KernelFunctor<cl::Buffer, cl::Buffer> * mul_kernel;
     cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer> *add_kernel;
     cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer> *triad_kernel;
+    cl::KernelFunctor<cl::Buffer, cl::Buffer, cl::Buffer, cl::LocalSpaceArg> *dot_kernel;
 
   public:
 
@@ -51,6 +59,7 @@ class OCLStream : public Stream<T>
     virtual void add() override;
     virtual void mul() override;
     virtual void triad() override;
+    virtual T dot() override;
 
     virtual void write_arrays(const std::vector<T>& a, const std::vector<T>& b, const std::vector<T>& c) override;
     virtual void read_arrays(std::vector<T>& a, std::vector<T>& b, std::vector<T>& c) override;