Added XGER routine, kernel, and tuner

6 files changed, 412 insertions, 8 deletions

diff --git a/src/clblast.cc b/src/clblast.cc
index 77999aaf..aed3f141 100644
--- a/src/clblast.cc
+++ b/src/clblast.cc
@@ -38,6 +38,7 @@
 #include "internal/routines/level2/xtrmv.h"
 #include "internal/routines/level2/xtbmv.h"
 #include "internal/routines/level2/xtpmv.h"
+#include "internal/routines/level2/xger.h"
 
 // BLAS level-3 includes
 #include "internal/routines/level3/xgemm.h"
@@ -835,14 +836,24 @@ template StatusCode Tpsv<double2>(const Layout, const Triangle, const Transpose,
 
 // General rank-1 matrix update: SGER/DGER
 template <typename T>
-StatusCode Ger(const Layout,
-               const size_t, const size_t,
-               const T,
-               const cl_mem, const size_t, const size_t,
-               const cl_mem, const size_t, const size_t,
-               cl_mem, const size_t, const size_t,
-               cl_command_queue*, cl_event*) {
-  return StatusCode::kNotImplemented;
+StatusCode Ger(const Layout layout,
+               const size_t m, const size_t n,
+               const T alpha,
+               const cl_mem x_buffer, const size_t x_offset, const size_t x_inc,
+               const cl_mem y_buffer, const size_t y_offset, const size_t y_inc,
+               cl_mem a_buffer, const size_t a_offset, const size_t a_ld,
+               cl_command_queue* queue, cl_event* event) {
+  auto queue_cpp = Queue(*queue);
+  auto event_cpp = Event(*event);
+  auto routine = Xger<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoGer(layout,
+                       m, n,
+                       alpha,
+                       Buffer<T>(x_buffer), x_offset, x_inc,
+                       Buffer<T>(y_buffer), y_offset, y_inc,
+                       Buffer<T>(a_buffer), a_offset, a_ld);
 }
 template StatusCode Ger<float>(const Layout,
                                const size_t, const size_t,
diff --git a/src/database.cc b/src/database.cc
index ba0a56d9..addd85d3 100644
--- a/src/database.cc
+++ b/src/database.cc
@@ -15,6 +15,7 @@
 #include "internal/database/xaxpy.h"
 #include "internal/database/xdot.h"
 #include "internal/database/xgemv.h"
+#include "internal/database/xger.h"
 #include "internal/database/xgemm.h"
 #include "internal/database/copy.h"
 #include "internal/database/pad.h"
@@ -31,6 +32,7 @@ const std::vector<Database::DatabaseEntry> Database::database = {
   XaxpySingle, XaxpyDouble, XaxpyComplexSingle, XaxpyComplexDouble,
   XdotSingle, XdotDouble, XdotComplexSingle, XdotComplexDouble,
   XgemvSingle, XgemvDouble, XgemvComplexSingle, XgemvComplexDouble,
+  XgerSingle, XgerDouble, XgerComplexSingle, XgerComplexDouble,
   XgemmSingle, XgemmDouble, XgemmComplexSingle, XgemmComplexDouble,
   CopySingle, CopyDouble, CopyComplexSingle, CopyComplexDouble,
   PadSingle, PadDouble, PadComplexSingle, PadComplexDouble,
diff --git a/src/kernels/common.opencl b/src/kernels/common.opencl
index f2a2e7a7..973c123e 100644
--- a/src/kernels/common.opencl
+++ b/src/kernels/common.opencl
@@ -147,6 +147,13 @@ R"(
   #define AXPBY(e, a, b, c, d) e = a*b + c*d
 #endif
 
+// The scalar GER function
+#if PRECISION == 3232 || PRECISION == 6464
+  #define GER(e, a, b, c, d) real ab; ab.x = MulReal(a,b); ab.y = MulImag(a,b); e.x = MulReal(ab,c) + d.x; e.y = MulImag(ab,c) + d.y
+#else
+  #define GER(e, a, b, c, d) e = a*b*c + d
+#endif
+
 // The complex conjugate operation for complex transforms
 #if PRECISION == 3232 || PRECISION == 6464
   #define COMPLEX_CONJUGATE(value) value.x = value.x; value.y = -value.y
diff --git a/src/kernels/level2/xger.opencl b/src/kernels/level2/xger.opencl
new file mode 100644
index 00000000..aa765b6c
--- /dev/null
+++ b/src/kernels/level2/xger.opencl
@@ -0,0 +1,149 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file contains the Xger kernel (generic version) for rank-1 matrix update.
+//
+// =================================================================================================
+
+// Enables loading of this file using the C++ pre-processor's #include (C++11 standard raw string
+// literal). Comment-out this line for syntax-highlighting when developing.
+R"(
+
+// =================================================================================================
+
+// Parameters set by the tuner or by the database. Here they are given a basic default value in case
+// this kernel file is used outside of the CLBlast library.
+
+#ifndef WGS1
+  #define WGS1 8    // The local work-group size in first dimension
+#endif
+#ifndef WGS2
+  #define WGS2 8    // The local work-group size in second dimension
+#endif
+#ifndef WPT
+  #define WPT 1     // The amount of work-per-thread in both dimensions
+#endif
+
+// =================================================================================================
+
+// Row-major version of the kernel
+__attribute__((reqd_work_group_size(WGS1, WGS2, 1)))
+__kernel void Xger(const int max_one, const int max_two, const real alpha,
+                   const __global real* restrict xgm, const int x_offset, const int x_inc,
+                   const __global real* ygm, const int y_offset, const int y_inc,
+                   __global real* restrict agm, const int a_offset, const int a_ld,
+                   const int is_rowmajor) {
+
+  // Register storage for X and Y
+  real xvalues[WPT];
+  real yvalues[WPT];
+
+  // Row-major version
+  if (is_rowmajor) {
+
+    // Loads the X-vector
+    #pragma unroll
+    for (int w=0; w<WPT; ++w) {
+      const int id2 = w*get_global_size(1) + get_global_id(1);
+      if (id2 < max_two) {
+        xvalues[w] = xgm[id2*x_inc + x_offset];
+      }
+    }
+
+    // Loads the Y-vector
+    #pragma unroll
+    for (int w=0; w<WPT; ++w) {
+      const int id1 = w*get_global_size(0) + get_global_id(0);
+      if (id1 < max_one) {
+        yvalues[w] = ygm[id1*y_inc + y_offset];
+      }
+    }
+
+    // Loops over the work per thread twice
+    #pragma unroll
+    for (int w1=0; w1<WPT; ++w1) {
+      #pragma unroll
+      for (int w2=0; w2<WPT; ++w2) {
+
+        // Global thread IDs
+        const int id1 = w1*get_global_size(0) + get_global_id(0);
+        const int id2 = w2*get_global_size(1) + get_global_id(1);
+
+        if (id1 < max_one && id2 < max_two) {
+
+          // Loads the current value of the A matrix
+          const int a_index = id2*a_ld + id1 + a_offset;
+          const real avalue = agm[a_index];
+
+          // Computes result = alpha * x[i] * y[j] + a[i][j]
+          real result;
+          GER(result, alpha, xvalues[w2], yvalues[w1], avalue);
+          
+          // Stores the final result
+          agm[a_index] = result;
+        }
+      }
+    }
+  }
+
+  // Col-major version
+  else {
+
+    // Loads the X-vector
+    #pragma unroll
+    for (int w=0; w<WPT; ++w) {
+      const int id1 = w*get_global_size(0) + get_global_id(0);
+      if (id1 < max_one) {
+        xvalues[w] = xgm[id1*x_inc + x_offset];
+      }
+    }
+
+    // Loads the Y-vector
+    #pragma unroll
+    for (int w=0; w<WPT; ++w) {
+      const int id2 = w*get_global_size(1) + get_global_id(1);
+      if (id2 < max_two) {
+        yvalues[w] = ygm[id2*y_inc + y_offset];
+      }
+    }
+
+    // Loops over the work per thread twice
+    #pragma unroll
+    for (int w1=0; w1<WPT; ++w1) {
+      #pragma unroll
+      for (int w2=0; w2<WPT; ++w2) {
+
+        // Global thread IDs
+        const int id1 = w1*get_global_size(0) + get_global_id(0);
+        const int id2 = w2*get_global_size(1) + get_global_id(1);
+
+        if (id1 < max_one && id2 < max_two) {
+
+          // Loads the current value of the A matrix
+          const int a_index = id2*a_ld + id1 + a_offset;
+          const real avalue = agm[a_index];
+
+          // Computes result = alpha * x[i] * y[j] + a[i][j]
+          real result;
+          GER(result, alpha, xvalues[w1], yvalues[w2], avalue);
+          
+          // Stores the final result
+          agm[a_index] = result;
+        }
+      }
+    }
+  }
+}
+
+// =================================================================================================
+
+// End of the C++11 raw string literal
+)"
+
+// =================================================================================================
diff --git a/src/routines/level2/xger.cc b/src/routines/level2/xger.cc
new file mode 100644
index 00000000..c3a24264
--- /dev/null
+++ b/src/routines/level2/xger.cc
@@ -0,0 +1,107 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file implements the Xger class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xger.h"
+
+#include <string>
+#include <vector>
+
+namespace clblast {
+// =================================================================================================
+
+// Specific implementations to get the memory-type based on a template argument
+template <> const Precision Xger<float>::precision_ = Precision::kSingle;
+template <> const Precision Xger<double>::precision_ = Precision::kDouble;
+
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xger<T>::Xger(Queue &queue, Event &event, const std::string &name):
+    Routine<T>(queue, event, name, {"Xger"}, precision_) {
+  source_string_ =
+    #include "../../kernels/level2/xger.opencl"
+  ;
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xger<T>::DoGer(const Layout layout,
+                          const size_t m, const size_t n,
+                          const T alpha,
+                          const Buffer<T> &x_buffer, const size_t x_offset, const size_t x_inc,
+                          const Buffer<T> &y_buffer, const size_t y_offset, const size_t y_inc,
+                          const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld) {
+
+  // Makes sure all dimensions are larger than zero
+  if (m == 0 || n == 0) { return StatusCode::kInvalidDimension; }
+
+  // Computes whether or not the matrix has an alternative layout (row or column-major).
+  const auto a_is_rowmajor = (layout == Layout::kRowMajor);
+  const auto a_one = (a_is_rowmajor) ? n : m;
+  const auto a_two = (a_is_rowmajor) ? m : n;
+
+  // Tests the matrix and the vectors for validity
+  auto status = TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld, sizeof(T));
+  if (ErrorIn(status)) { return status; }
+  status = TestVectorX(m, x_buffer, x_offset, x_inc, sizeof(T));
+  if (ErrorIn(status)) { return status; }
+  status = TestVectorY(n, y_buffer, y_offset, y_inc, sizeof(T));
+  if (ErrorIn(status)) { return status; }
+
+  // Retrieves the Xgemv kernel from the compiled binary
+  try {
+    auto& program = GetProgramFromCache();
+    auto kernel = Kernel(program, "Xger");
+
+    // Sets the kernel arguments
+    kernel.SetArgument(0, static_cast<int>(a_one));
+    kernel.SetArgument(1, static_cast<int>(a_two));
+    kernel.SetArgument(2, alpha);
+    kernel.SetArgument(3, x_buffer());
+    kernel.SetArgument(4, static_cast<int>(x_offset));
+    kernel.SetArgument(5, static_cast<int>(x_inc));
+    kernel.SetArgument(6, y_buffer());
+    kernel.SetArgument(7, static_cast<int>(y_offset));
+    kernel.SetArgument(8, static_cast<int>(y_inc));
+    kernel.SetArgument(9, a_buffer());
+    kernel.SetArgument(10, static_cast<int>(a_offset));
+    kernel.SetArgument(11, static_cast<int>(a_ld));
+    kernel.SetArgument(12, static_cast<int>(a_is_rowmajor));
+
+    // Launches the kernel
+    auto a_one_ceiled = CeilDiv(Ceil(a_one, db_["WGS1"]), db_["WPT"]);
+    auto a_two_ceiled = CeilDiv(Ceil(a_two, db_["WGS2"]), db_["WPT"]);
+    auto global = std::vector<size_t>{a_one_ceiled, a_two_ceiled};
+    auto local = std::vector<size_t>{db_["WGS1"], db_["WGS2"]};
+    status = RunKernel(kernel, global, local);
+    if (ErrorIn(status)) { return status; }
+
+    // Waits for all kernels to finish
+    queue_.Finish();
+
+    // Succesfully finished the computation
+    return StatusCode::kSuccess;
+  } catch (...) { return StatusCode::kInvalidKernel; }
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xger<float>;
+template class Xger<double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/tuning/xger.cc b/src/tuning/xger.cc
new file mode 100644
index 00000000..a47f5a78
--- /dev/null
+++ b/src/tuning/xger.cc
@@ -0,0 +1,128 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file uses the CLTune auto-tuner to tune the xger OpenCL kernels.
+//
+// =================================================================================================
+
+#include <string>
+#include <vector>
+
+#include "internal/utilities.h"
+#include "internal/tuning.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class TuneXger {
+ public:
+
+  // The representative kernel and the source code
+  static std::string KernelFamily() { return "xger"; }
+  static std::string KernelName() { return "Xger"; }
+  static std::string GetSources() {
+    return
+      #include "../src/kernels/common.opencl"
+      #include "../src/kernels/level2/xger.opencl"
+    ;
+  }
+
+  // The list of arguments relevant for this routine
+  static std::vector<std::string> GetOptions() { return {kArgN, kArgM, kArgAlpha}; }
+
+  // Tests for valid arguments
+  static void TestValidArguments(const Arguments<T> &) { }
+
+  // Sets the default values for the arguments
+  static size_t DefaultM() { return 1024; }
+  static size_t DefaultN() { return 1024; }
+  static size_t DefaultK() { return 1; } // N/A for this kernel
+  static double DefaultFraction() { return 1.0; } // N/A for this kernel
+
+  // Describes how to obtain the sizes of the buffers
+  static size_t GetSizeX(const Arguments<T> &args) { return args.m; }
+  static size_t GetSizeY(const Arguments<T> &args) { return args.n; }
+  static size_t GetSizeA(const Arguments<T> &args) { return args.m * args.n; }
+  static size_t GetSizeB(const Arguments<T> &) { return 1; } // N/A for this kernel
+  static size_t GetSizeC(const Arguments<T> &) { return 1; } // N/A for this kernel
+  static size_t GetSizeTemp(const Arguments<T> &) { return 1; } // N/A for this kernel
+
+  // Sets the tuning parameters and their possible values
+  static void SetParameters(cltune::Tuner &tuner, const size_t id) {
+    tuner.AddParameter(id, "WGS1", {4, 8, 16, 32, 64, 128, 256, 512});
+    tuner.AddParameter(id, "WGS2", {1, 2, 4, 8, 16, 32, 64, 128, 256});
+    tuner.AddParameter(id, "WPT", {1, 2, 4});
+  }
+
+  // Sets the constraints and local memory size
+  static void SetConstraints(cltune::Tuner &, const size_t) { }
+  static void SetLocalMemorySize(cltune::Tuner &, const size_t, const Arguments<T> &) { }
+
+  // Sets the base thread configuration
+  static std::vector<size_t> GlobalSize(const Arguments<T> &args) { return {args.m, args.n}; }
+  static std::vector<size_t> GlobalSizeRef(const Arguments<T> &args) { return GlobalSize(args); }
+  static std::vector<size_t> LocalSize() { return {1, 1}; }
+  static std::vector<size_t> LocalSizeRef() { return {8, 8}; }
+
+  // Transforms the thread configuration based on the parameters
+  using TransformVector = std::vector<std::vector<std::string>>;
+  static TransformVector MulLocal() { return {{"WGS1", "WGS2"}}; }
+  static TransformVector DivLocal() { return {}; }
+  static TransformVector MulGlobal() { return {}; }
+  static TransformVector DivGlobal() { return {{"WPT", "WPT"}}; }
+
+  // Sets the kernel's arguments
+  static void SetArguments(cltune::Tuner &tuner, const Arguments<T> &args,
+                           std::vector<T> &x_vec, std::vector<T> &y_vec,
+                           std::vector<T> &a_mat, std::vector<T> &, std::vector<T> &,
+                           std::vector<T> &) {
+    tuner.AddArgumentScalar(static_cast<int>(args.m));
+    tuner.AddArgumentScalar(static_cast<int>(args.n));
+    tuner.AddArgumentScalar(args.alpha);
+    tuner.AddArgumentInput(x_vec);
+    tuner.AddArgumentScalar(0); // x_offset
+    tuner.AddArgumentScalar(1); // x_increment
+    tuner.AddArgumentInput(y_vec);
+    tuner.AddArgumentScalar(0); // y_offset
+    tuner.AddArgumentScalar(1); // y_increment
+    tuner.AddArgumentOutput(a_mat);
+    tuner.AddArgumentScalar(0); // a_offset
+    tuner.AddArgumentScalar(static_cast<int>(args.m)); // a_ld
+    tuner.AddArgumentScalar(0); // a_is_rowmajor
+  }
+
+  // Describes how to compute the performance metrics
+  static size_t GetMetric(const Arguments<T> &args) {
+    return (2*args.m*args.n + args.m + args.n) * GetBytes(args.precision);
+  }
+  static std::string PerformanceUnit() { return "GB/s"; }
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// Shortcuts to the clblast namespace
+using float2 = clblast::float2;
+using double2 = clblast::double2;
+
+// Main function (not within the clblast namespace)
+int main(int argc, char *argv[]) {
+  switch(clblast::GetPrecision(argc, argv)) {
+    case clblast::Precision::kHalf: throw std::runtime_error("Unsupported precision mode");
+    case clblast::Precision::kSingle: clblast::Tuner<clblast::TuneXger<float>, float>(argc, argv); break;
+    case clblast::Precision::kDouble: clblast::Tuner<clblast::TuneXger<double>, double>(argc, argv); break;
+    case clblast::Precision::kComplexSingle: clblast::Tuner<clblast::TuneXger<float2>, float2>(argc, argv); break;
+    case clblast::Precision::kComplexDouble: clblast::Tuner<clblast::TuneXger<double2>, double2>(argc, argv); break;
+  }
+  return 0;
+}
+
+// =================================================================================================