Added support for xHER, xHPR, xSYR, and xSPR routines

10 files changed, 519 insertions, 90 deletions

diff --git a/src/clblast.cc b/src/clblast.cc
index 47f2c59d..466de83e 100644
--- a/src/clblast.cc
+++ b/src/clblast.cc
@@ -41,6 +41,10 @@
 #include "internal/routines/level2/xger.h"
 #include "internal/routines/level2/xgeru.h"
 #include "internal/routines/level2/xgerc.h"
+#include "internal/routines/level2/xher.h"
+#include "internal/routines/level2/xhpr.h"
+#include "internal/routines/level2/xsyr.h"
+#include "internal/routines/level2/xspr.h"
 
 // BLAS level-3 includes
 #include "internal/routines/level3/xgemm.h"
@@ -946,13 +950,22 @@ template StatusCode Gerc<double2>(const Layout,
 
 // Hermitian rank-1 matrix update: CHER/ZHER
 template <typename T>
-StatusCode Her(const Layout, const Triangle,
-               const size_t,
-               const 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 Her(const Layout layout, const Triangle triangle,
+               const size_t n,
+               const T alpha,
+               const cl_mem x_buffer, const size_t x_offset, const size_t x_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 = Xher<std::complex<T>,T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoHer(layout, triangle,
+                       n,
+                       alpha,
+                       Buffer<std::complex<T>>(x_buffer), x_offset, x_inc,
+                       Buffer<std::complex<T>>(a_buffer), a_offset, a_ld);
 }
 template StatusCode Her<float>(const Layout, const Triangle,
                                const size_t,
@@ -969,13 +982,22 @@ template StatusCode Her<double>(const Layout, const Triangle,
 
 // Hermitian packed rank-1 matrix update: CHPR/ZHPR
 template <typename T>
-StatusCode Hpr(const Layout, const Triangle,
-               const size_t,
-               const T,
-               const cl_mem, const size_t, const size_t,
-               cl_mem, const size_t,
-               cl_command_queue*, cl_event*) {
-  return StatusCode::kNotImplemented;
+StatusCode Hpr(const Layout layout, const Triangle triangle,
+               const size_t n,
+               const T alpha,
+               const cl_mem x_buffer, const size_t x_offset, const size_t x_inc,
+               cl_mem ap_buffer, const size_t ap_offset,
+               cl_command_queue* queue, cl_event* event) {
+  auto queue_cpp = Queue(*queue);
+  auto event_cpp = Event(*event);
+  auto routine = Xhpr<std::complex<T>,T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoHpr(layout, triangle,
+                       n,
+                       alpha,
+                       Buffer<std::complex<T>>(x_buffer), x_offset, x_inc,
+                       Buffer<std::complex<T>>(ap_buffer), ap_offset);
 }
 template StatusCode Hpr<float>(const Layout, const Triangle,
                                const size_t,
@@ -1044,13 +1066,22 @@ template StatusCode Hpr2<double2>(const Layout, const Triangle,
 
 // Symmetric rank-1 matrix update: SSYR/DSYR
 template <typename T>
-StatusCode Syr(const Layout, const Triangle,
-               const size_t,
-               const 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 Syr(const Layout layout, const Triangle triangle,
+               const size_t n,
+               const T alpha,
+               const cl_mem x_buffer, const size_t x_offset, const size_t x_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 = Xsyr<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoSyr(layout, triangle,
+                       n,
+                       alpha,
+                       Buffer<T>(x_buffer), x_offset, x_inc,
+                       Buffer<T>(a_buffer), a_offset, a_ld);
 }
 template StatusCode Syr<float>(const Layout, const Triangle,
                                const size_t,
@@ -1067,13 +1098,22 @@ template StatusCode Syr<double>(const Layout, const Triangle,
 
 // Symmetric packed rank-1 matrix update: SSPR/DSPR
 template <typename T>
-StatusCode Spr(const Layout, const Triangle,
-               const size_t,
-               const T,
-               const cl_mem, const size_t, const size_t,
-               cl_mem, const size_t,
-               cl_command_queue*, cl_event*) {
-  return StatusCode::kNotImplemented;
+StatusCode Spr(const Layout layout, const Triangle triangle,
+               const size_t n,
+               const T alpha,
+               const cl_mem x_buffer, const size_t x_offset, const size_t x_inc,
+               cl_mem ap_buffer, const size_t ap_offset,
+               cl_command_queue* queue, cl_event* event) {
+  auto queue_cpp = Queue(*queue);
+  auto event_cpp = Event(*event);
+  auto routine = Xspr<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoSpr(layout, triangle,
+                       n,
+                       alpha,
+                       Buffer<T>(x_buffer), x_offset, x_inc,
+                       Buffer<T>(ap_buffer), ap_offset);
 }
 template StatusCode Spr<float>(const Layout, const Triangle,
                                const size_t,
diff --git a/src/kernels/level2/level2.opencl b/src/kernels/level2/level2.opencl
new file mode 100644
index 00000000..ea7d7afd
--- /dev/null
+++ b/src/kernels/level2/level2.opencl
@@ -0,0 +1,94 @@
+
+// =================================================================================================
+// 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 common functions for matrix update kernels (Xger, Xher).
+//
+// =================================================================================================
+
+// 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
+
+// =================================================================================================
+
+// Returns an element from a vector
+inline real LoadVector(const int id, const int max,
+                       __global real* restrict gm, const int offset, const int inc,
+                       const int do_conjugate) {
+  if (id < max) {
+    real result = gm[id*inc + offset];
+    if (do_conjugate) {
+      #if defined(ROUTINE_GERC)
+        COMPLEX_CONJUGATE(result);
+      #endif
+      #if defined(ROUTINE_HER)
+        COMPLEX_CONJUGATE(result);
+      #endif
+    }
+    return result;
+  }
+  else {
+    real default_result;
+    SetToZero(default_result);
+    return default_result;
+  }
+}
+
+// Performs the rank-1 matrix update
+inline void MatrixUpdate(const int id1, const int id2, const int max1, const int max2,
+                         __global real* restrict agm, const int a_offset, const int a_ld,
+                         const real alpha, const real xvalue, const real yvalue) {
+
+  // Bounds of a regular matrix
+  if (id1 < max1 && id2 < max2) {
+
+    #if defined(ROUTINE_SPR) || defined(ROUTINE_HPR)
+      const int a_index = (id1 <= id2) ? ((id2+1)*id2)/2 + id1 + a_offset : ((id1+1)*id1)/2 + id2 + a_offset;
+    #else
+      const int a_index = id2*a_ld + id1 + a_offset;
+    #endif
+
+    // Loads the current value of the A matrix
+    const real avalue = agm[a_index];
+
+    // Computes result = alpha * x[i] * y[j] + a[i][j]
+    real result;
+    GER(result, alpha, xvalue, yvalue, avalue);
+
+    // For hermetian matrices
+    #if defined(ROUTINE_HER) || defined(ROUTINE_HPR)
+      if (id1 == id2) { result.y = ZERO; }
+    #endif
+    
+    // Stores the final result
+    agm[a_index] = result;
+  }
+}
+
+// =================================================================================================
+
+// End of the C++11 raw string literal
+)"
+
+// =================================================================================================
diff --git a/src/kernels/level2/xger.opencl b/src/kernels/level2/xger.opencl
index c1cddea1..ce8e04bb 100644
--- a/src/kernels/level2/xger.opencl
+++ b/src/kernels/level2/xger.opencl
@@ -7,7 +7,7 @@
 // Author(s):
 //   Cedric Nugteren <www.cedricnugteren.nl>
 //
-// This file contains the Xger kernel (generic version) for rank-1 matrix update.
+// This file contains the Xger kernels for rank-1 matrix update.
 //
 // =================================================================================================
 
@@ -17,24 +17,9 @@ 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
+// Regular version of the rank-1 matrix update kernel (GER, GERU, GERC)
 __attribute__((reqd_work_group_size(WGS1, WGS2, 1)))
-__kernel void Xger(const int max_one, const int max_two, const real alpha,
+__kernel void Xger(const int max1, const int max2, 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,
@@ -51,21 +36,14 @@ __kernel void Xger(const int max_one, const int max_two, const real alpha,
     #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];
-      }
+      xvalues[w] = LoadVector(id2, max2, xgm, x_offset, x_inc, false);
     }
 
     // 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];
-        #if defined(ROUTINE_GERC)
-          COMPLEX_CONJUGATE(yvalues[w]);
-        #endif
-      }
+      yvalues[w] = LoadVector(id1, max1, ygm, y_offset, y_inc, true);
     }
 
     // Loops over the work per thread twice
@@ -78,19 +56,9 @@ __kernel void Xger(const int max_one, const int max_two, const real alpha,
         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;
-        }
+        // Loads A, performs the operation, and stores the result into A
+        MatrixUpdate(id1, id2, max1, max2, agm, a_offset, a_ld,
+                     alpha, xvalues[w2], yvalues[w1]);
       }
     }
   }
@@ -102,21 +70,14 @@ __kernel void Xger(const int max_one, const int max_two, const real alpha,
     #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];
-      }
+      xvalues[w] = LoadVector(id1, max1, xgm, x_offset, x_inc, false);
     }
 
     // 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];
-        #if defined(ROUTINE_GERC)
-          COMPLEX_CONJUGATE(yvalues[w]);
-        #endif
-      }
+      yvalues[w] = LoadVector(id2, max2, ygm, y_offset, y_inc, true);
     }
 
     // Loops over the work per thread twice
@@ -129,19 +90,9 @@ __kernel void Xger(const int max_one, const int max_two, const real alpha,
         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;
-        }
+        // Loads A, performs the operation, and stores the result into A
+        MatrixUpdate(id1, id2, max1, max2, agm, a_offset, a_ld,
+                     alpha, xvalues[w1], yvalues[w2]);
       }
     }
   }
diff --git a/src/kernels/level2/xher.opencl b/src/kernels/level2/xher.opencl
new file mode 100644
index 00000000..13bc4135
--- /dev/null
+++ b/src/kernels/level2/xher.opencl
@@ -0,0 +1,73 @@
+
+// =================================================================================================
+// 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 Xher kernels 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"(
+
+// =================================================================================================
+
+// Symmetric version of the rank-1 matrix update kernel (HER, HPR, SYR, SPR)
+__attribute__((reqd_work_group_size(WGS1, WGS2, 1)))
+__kernel void Xher(const int n, const real alpha,
+                   const __global real* restrict xgm, const int x_offset, const int x_inc,
+                   __global real* restrict agm, const int a_offset, const int a_ld,
+                   const int is_upper, const int is_rowmajor) {
+
+  // Register storage for X and XT
+  real xvalues[WPT];
+  real xtvalues[WPT];
+
+  // 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);
+    xvalues[w] = LoadVector(id2, n, xgm, x_offset, x_inc, !is_rowmajor);
+  }
+
+  // Loads the X-transposed-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id1 = w*get_global_size(0) + get_global_id(0);
+    xtvalues[w] = LoadVector(id1, n, xgm, x_offset, x_inc, is_rowmajor);
+  }
+
+  // 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);
+
+      // Skip these threads if they do not contain threads contributing to the matrix-triangle
+      if ((is_upper && (id1 > id2)) || (!is_upper && (id2 > id1))) {
+        // Do nothing
+      }
+
+      // Loads A, performs the operation, and stores the result into A
+      else {
+        MatrixUpdate(id1, id2, n, n, agm, a_offset, a_ld, alpha, xvalues[w2], xtvalues[w1]);
+      }
+    }
+  }
+}
+
+// =================================================================================================
+
+// End of the C++11 raw string literal
+)"
+
+// =================================================================================================
diff --git a/src/routines/level2/xger.cc b/src/routines/level2/xger.cc
index 55fa26d4..0953c8bb 100644
--- a/src/routines/level2/xger.cc
+++ b/src/routines/level2/xger.cc
@@ -32,6 +32,7 @@ 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/level2.opencl"
     #include "../../kernels/level2/xger.opencl"
   ;
 }
diff --git a/src/routines/level2/xher.cc b/src/routines/level2/xher.cc
new file mode 100644
index 00000000..4227252e
--- /dev/null
+++ b/src/routines/level2/xher.cc
@@ -0,0 +1,119 @@
+
+// =================================================================================================
+// 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 Xher class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xher.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Specific implementations to get the memory-type based on a template argument
+template <> const Precision Xher<float, float>::precision_ = Precision::kSingle;
+template <> const Precision Xher<double, double>::precision_ = Precision::kDouble;
+template <> const Precision Xher<float2, float>::precision_ = Precision::kComplexSingle;
+template <> const Precision Xher<double2, double>::precision_ = Precision::kComplexDouble;
+
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T, typename U>
+Xher<T,U>::Xher(Queue &queue, Event &event, const std::string &name):
+    Routine<T>(queue, event, name, {"Xger"}, precision_) {
+  source_string_ =
+    #include "../../kernels/level2/level2.opencl"
+    #include "../../kernels/level2/xher.opencl"
+  ;
+}
+
+// =================================================================================================
+
+// Specializations to compute alpha of type 'T'
+template <> float2 Xher<float2,float>::GetAlpha(const float alpha) { return float2{alpha, 0.0f}; }
+template <> double2 Xher<double2,double>::GetAlpha(const double alpha) { return double2{alpha, 0.0}; }
+template <> float Xher<float,float>::GetAlpha(const float alpha) { return alpha; }
+template <> double Xher<double,double>::GetAlpha(const double alpha) { return alpha; }
+
+// =================================================================================================
+
+// The main routine
+template <typename T, typename U>
+StatusCode Xher<T,U>::DoHer(const Layout layout, const Triangle triangle,
+                            const size_t n,
+                            const U alpha,
+                            const Buffer<T> &x_buffer, const size_t x_offset, const size_t x_inc,
+                            const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+                            const bool packed) {
+
+  // Makes sure the dimensions are larger than zero
+  if (n == 0) { return StatusCode::kInvalidDimension; }
+
+  // The data is either in the upper or lower triangle
+  const auto is_upper = ((triangle == Triangle::kUpper && layout != Layout::kRowMajor) ||
+                         (triangle == Triangle::kLower && layout == Layout::kRowMajor));
+  const auto is_rowmajor = (layout == Layout::kRowMajor);
+
+  // Creates a matching version of alpha
+  const auto matching_alpha = GetAlpha(alpha);
+
+  // Tests the matrix and the vectors for validity
+  auto status = StatusCode::kSuccess;
+  if (packed) { status = TestMatrixAP(n, a_buffer, a_offset, sizeof(T)); }
+  else { status = TestMatrixA(n, n, a_buffer, a_offset, a_ld, sizeof(T)); }
+  if (ErrorIn(status)) { return status; }
+  status = TestVectorX(n, x_buffer, x_offset, x_inc, sizeof(T));
+  if (ErrorIn(status)) { return status; }
+
+  // Retrieves the Xgemv kernel from the compiled binary
+  try {
+    auto& program = GetProgramFromCache();
+    auto kernel = Kernel(program, "Xher");
+
+    // Sets the kernel arguments
+    kernel.SetArgument(0, static_cast<int>(n));
+    kernel.SetArgument(1, matching_alpha);
+    kernel.SetArgument(2, x_buffer());
+    kernel.SetArgument(3, static_cast<int>(x_offset));
+    kernel.SetArgument(4, static_cast<int>(x_inc));
+    kernel.SetArgument(5, a_buffer());
+    kernel.SetArgument(6, static_cast<int>(a_offset));
+    kernel.SetArgument(7, static_cast<int>(a_ld));
+    kernel.SetArgument(8, static_cast<int>(is_upper));
+    kernel.SetArgument(9, static_cast<int>(is_rowmajor));
+
+    // Launches the kernel
+    auto global_one = CeilDiv(Ceil(n, db_["WGS1"]), db_["WPT"]);
+    auto global_two = CeilDiv(Ceil(n, db_["WGS2"]), db_["WPT"]);
+    auto global = std::vector<size_t>{global_one, global_two};
+    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 Xher<float, float>;
+template class Xher<double, double>;
+template class Xher<float2, float>;
+template class Xher<double2, double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xhpr.cc b/src/routines/level2/xhpr.cc
new file mode 100644
index 00000000..abe00669
--- /dev/null
+++ b/src/routines/level2/xhpr.cc
@@ -0,0 +1,50 @@
+
+// =================================================================================================
+// 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 Xhpr class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xhpr.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T, typename U>
+Xhpr<T,U>::Xhpr(Queue &queue, Event &event, const std::string &name):
+    Xher<T,U>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T, typename U>
+StatusCode Xhpr<T,U>::DoHpr(const Layout layout, const Triangle triangle,
+                            const size_t n,
+                            const U alpha,
+                            const Buffer<T> &x_buffer, const size_t x_offset, const size_t x_inc,
+                            const Buffer<T> &ap_buffer, const size_t ap_offset) {
+
+  // 
+  return DoHer(layout, triangle, n, alpha,
+               x_buffer, x_offset, x_inc,
+               ap_buffer, ap_offset, n, true);
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xhpr<float2, float>;
+template class Xhpr<double2, double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xspr.cc b/src/routines/level2/xspr.cc
new file mode 100644
index 00000000..5159ad50
--- /dev/null
+++ b/src/routines/level2/xspr.cc
@@ -0,0 +1,50 @@
+
+// =================================================================================================
+// 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 Xspr class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xspr.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xspr<T>::Xspr(Queue &queue, Event &event, const std::string &name):
+    Xher<T,T>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xspr<T>::DoSpr(const Layout layout, const Triangle triangle,
+                          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> &ap_buffer, const size_t ap_offset) {
+
+  // 
+  return DoHer(layout, triangle, n, alpha,
+               x_buffer, x_offset, x_inc,
+               ap_buffer, ap_offset, n, true);
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xspr<float>;
+template class Xspr<double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xsyr.cc b/src/routines/level2/xsyr.cc
new file mode 100644
index 00000000..755fde0d
--- /dev/null
+++ b/src/routines/level2/xsyr.cc
@@ -0,0 +1,50 @@
+
+// =================================================================================================
+// 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 Xsyr class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xsyr.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xsyr<T>::Xsyr(Queue &queue, Event &event, const std::string &name):
+    Xher<T,T>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xsyr<T>::DoSyr(const Layout layout, const Triangle triangle,
+                          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> &a_buffer, const size_t a_offset, const size_t a_ld) {
+
+  // 
+  return DoHer(layout, triangle, n, alpha,
+               x_buffer, x_offset, x_inc,
+               a_buffer, a_offset, a_ld);
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xsyr<float>;
+template class Xsyr<double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/tuning/xger.cc b/src/tuning/xger.cc
index a47f5a78..39efdb81 100644
--- a/src/tuning/xger.cc
+++ b/src/tuning/xger.cc
@@ -31,6 +31,7 @@ class TuneXger {
   static std::string GetSources() {
     return
       #include "../src/kernels/common.opencl"
+      #include "../src/kernels/level2/level2.opencl"
       #include "../src/kernels/level2/xger.opencl"
     ;
   }