Merge branch 'rank2_update_routines' into development

23 files changed, 1245 insertions, 64 deletions

diff --git a/CHANGELOG b/CHANGELOG
index 72f53550..3c91c31b 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -9,8 +9,12 @@ Development version (next release)
   * CGERC/ZGERC
   * CHER/ZHER
   * CHPR/ZHPR
+  * CHER2/ZHER2
+  * CHPR2/ZHPR2
   * CSYR/ZSYR
   * CSPR/ZSPR
+  * CSYR2/ZSYR2
+  * CSPR2/ZSPR2
 
 Version 0.5.0
 - Improved structure and performance of level-2 routines (xSYMV/xHEMV)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e35fda7a..508dad72 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -112,7 +112,7 @@ set(SAMPLE_PROGRAMS_CPP sgemm)
 set(SAMPLE_PROGRAMS_C sgemm)
 set(LEVEL1_ROUTINES xswap xscal xcopy xaxpy xdot xdotu xdotc)
 set(LEVEL2_ROUTINES xgemv xgbmv xhemv xhbmv xhpmv xsymv xsbmv xspmv xtrmv xtbmv xtpmv
-                    xger xgeru xgerc xher xhpr xsyr xspr)
+                    xger xgeru xgerc xher xhpr xher2 xhpr2 xsyr xspr xsyr2 xspr2)
 set(LEVEL3_ROUTINES xgemm xsymm xhemm xsyrk xherk xsyr2k xher2k xtrmm)
 set(ROUTINES ${LEVEL1_ROUTINES} ${LEVEL2_ROUTINES} ${LEVEL3_ROUTINES})
 set(PRECISIONS 32 64 3232 6464)
diff --git a/README.md b/README.md
index 3da0e262..87c7bb53 100644
--- a/README.md
+++ b/README.md
@@ -176,7 +176,6 @@ CLBlast is in active development but already supports almost all the BLAS routin
 | xASUM    |   |   | - | - | +SC +DZ |
 | IxAMAX   |   |   |   |   |         |
 
-
 | Level-2  | S | D | C | Z | Notes   |
 | ---------|---|---|---|---|---------|
 | xGEMV    | ✔ | ✔ | ✔ | ✔ |         |
@@ -198,12 +197,12 @@ CLBlast is in active development but already supports almost all the BLAS routin
 | xGERC    | - | - | ✔ | ✔ |         |
 | xHER     | - | - | ✔ | ✔ |         |
 | xHPR     | - | - | ✔ | ✔ |         |
-| xHER2    | - | - |   |   |         |
-| xHPR2    | - | - |   |   |         |
+| xHER2    | - | - | ✔ | ✔ |         |
+| xHPR2    | - | - | ✔ | ✔ |         |
 | xSYR     | ✔ | ✔ | - | - |         |
 | xSPR     | ✔ | ✔ | - | - |         |
-| xSYR2    |   |   | - | - |         |
-| xSPR2    |   |   | - | - |         |
+| xSYR2    | ✔ | ✔ | - | - |         |
+| xSPR2    | ✔ | ✔ | - | - |         |
 
 | Level-3  | S | D | C | Z | Notes   |
 | ---------|---|---|---|---|---------|
@@ -246,11 +245,8 @@ To-do list before release of version 1.0
 - Improve device performance:
   * Tune for a wider range of devices
   * Allow users to define custom tuned parameters
-- Improve the tuning
-  * Make the tuners upload their data to a central server
 - Improve the performance comparisons:
   * Enable comparison against optionally: ViennaCL, cuBLAS, MAGMA OpenCL
 - Further reduce the likelihood of crashes:
   * Add checks for proper command-line arguments in the tuner, tester and client
-  * Add checks for valid database parameters
   * Test in multi-threaded environments
diff --git a/include/internal/routines/level2/xher2.h b/include/internal/routines/level2/xher2.h
new file mode 100644
index 00000000..26f69046
--- /dev/null
+++ b/include/internal/routines/level2/xher2.h
@@ -0,0 +1,60 @@
+
+// =================================================================================================
+// 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 Xher2 routine. The precision is implemented using a template argument.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_ROUTINES_XHER2_H_
+#define CLBLAST_ROUTINES_XHER2_H_
+
+#include "internal/routine.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class Xher2: public Routine<T> {
+ public:
+
+  // Members and methods from the base class
+  using Routine<T>::db_;
+  using Routine<T>::source_string_;
+  using Routine<T>::queue_;
+  using Routine<T>::GetProgramFromCache;
+  using Routine<T>::TestVectorX;
+  using Routine<T>::TestVectorY;
+  using Routine<T>::TestMatrixA;
+  using Routine<T>::TestMatrixAP;
+  using Routine<T>::RunKernel;
+  using Routine<T>::ErrorIn;
+
+  // Constructor
+  Xher2(Queue &queue, Event &event, const std::string &name = "HER2");
+
+  // Templated-precision implementation of the routine
+  StatusCode DoHer2(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> &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,
+                    const bool packed = false);
+
+ private:
+  // Static variable to get the precision
+  const static Precision precision_;
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_ROUTINES_XHER2_H_
+#endif
diff --git a/include/internal/routines/level2/xhpr2.h b/include/internal/routines/level2/xhpr2.h
new file mode 100644
index 00000000..fd243d33
--- /dev/null
+++ b/include/internal/routines/level2/xhpr2.h
@@ -0,0 +1,46 @@
+
+// =================================================================================================
+// 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 Xhpr2 routine. The precision is implemented using a template argument.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_ROUTINES_XHPR2_H_
+#define CLBLAST_ROUTINES_XHPR2_H_
+
+#include "internal/routines/level2/xher2.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class Xhpr2: public Xher2<T> {
+ public:
+
+  // Uses the regular Xher2 routine
+  using Xher2<T>::DoHer2;
+
+  // Constructor
+  Xhpr2(Queue &queue, Event &event, const std::string &name = "HPR2");
+
+  // Templated-precision implementation of the routine
+  StatusCode DoHpr2(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> &y_buffer, const size_t y_offset, const size_t y_inc,
+                    const Buffer<T> &ap_buffer, const size_t ap_offset);
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_ROUTINES_XHPR2_H_
+#endif
diff --git a/include/internal/routines/level2/xspr2.h b/include/internal/routines/level2/xspr2.h
new file mode 100644
index 00000000..3d5f4992
--- /dev/null
+++ b/include/internal/routines/level2/xspr2.h
@@ -0,0 +1,46 @@
+
+// =================================================================================================
+// 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 Xspr2 routine. The precision is implemented using a template argument.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_ROUTINES_XSPR2_H_
+#define CLBLAST_ROUTINES_XSPR2_H_
+
+#include "internal/routines/level2/xher2.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class Xspr2: public Xher2<T> {
+ public:
+
+  // Uses the regular Xher2 routine
+  using Xher2<T>::DoHer2;
+
+  // Constructor
+  Xspr2(Queue &queue, Event &event, const std::string &name = "SPR2");
+
+  // Templated-precision implementation of the routine
+  StatusCode DoSpr2(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> &y_buffer, const size_t y_offset, const size_t y_inc,
+                    const Buffer<T> &ap_buffer, const size_t ap_offset);
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_ROUTINES_XSPR2_H_
+#endif
diff --git a/include/internal/routines/level2/xsyr2.h b/include/internal/routines/level2/xsyr2.h
new file mode 100644
index 00000000..f4dc9375
--- /dev/null
+++ b/include/internal/routines/level2/xsyr2.h
@@ -0,0 +1,46 @@
+
+// =================================================================================================
+// 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 Xsyr2 routine. The precision is implemented using a template argument.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_ROUTINES_XSYR2_H_
+#define CLBLAST_ROUTINES_XSYR2_H_
+
+#include "internal/routines/level2/xher2.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class Xsyr2: public Xher2<T> {
+ public:
+
+  // Uses the regular Xher2 routine
+  using Xher2<T>::DoHer2;
+
+  // Constructor
+  Xsyr2(Queue &queue, Event &event, const std::string &name = "SYR2");
+
+  // Templated-precision implementation of the routine
+  StatusCode DoSyr2(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> &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);
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_ROUTINES_XSYR2_H_
+#endif
diff --git a/scripts/generator/generator.py b/scripts/generator/generator.py
index 0f5fbfa7..99f326cd 100644
--- a/scripts/generator/generator.py
+++ b/scripts/generator/generator.py
@@ -83,12 +83,12 @@ routines = [
   Routine(True,  "2b", "gerc",  T,  [C,Z],     ["m","n"], ["layout"], ["x","y"], ["a"], ["alpha"], False, "General rank-1 complex conjugated matrix update"),
   Routine(True,  "2b", "her",   Tc, [Css,Zdd], ["n"], ["layout","triangle"], ["x"], ["a"], ["alpha"], False, "Hermitian rank-1 matrix update"),
   Routine(True,  "2b", "hpr",   Tc, [Css,Zdd], ["n"], ["layout","triangle"], ["x"], ["ap"], ["alpha"], False, "Hermitian packed rank-1 matrix update"),
-  Routine(False, "2b", "her2",  T,  [C,Z],     ["n"], ["layout","triangle"], ["x","y"], ["a"], ["alpha"], False, "Hermitian rank-2 matrix update"),
-  Routine(False, "2b", "hpr2",  T,  [C,Z],     ["n"], ["layout","triangle"], ["x","y"], ["ap"], ["alpha"], False, "Hermitian packed rank-2 matrix update"),
+  Routine(True,  "2b", "her2",  T,  [C,Z],     ["n"], ["layout","triangle"], ["x","y"], ["a"], ["alpha"], False, "Hermitian rank-2 matrix update"),
+  Routine(True,  "2b", "hpr2",  T,  [C,Z],     ["n"], ["layout","triangle"], ["x","y"], ["ap"], ["alpha"], False, "Hermitian packed rank-2 matrix update"),
   Routine(True,  "2b", "syr",   T,  [S,D],     ["n"], ["layout","triangle"], ["x"], ["a"], ["alpha"], False, "Symmetric rank-1 matrix update"),
   Routine(True,  "2b", "spr",   T,  [S,D],     ["n"], ["layout","triangle"], ["x"], ["ap"], ["alpha"], False, "Symmetric packed rank-1 matrix update"),
-  Routine(False, "2b", "syr2",  T,  [S,D],     ["n"], ["layout","triangle"], ["x","y"], ["a"], ["alpha"], False, "Symmetric rank-2 matrix update"),
-  Routine(False, "2b", "spr2",  T,  [S,D],     ["n"], ["layout","triangle"], ["x","y"], ["ap"], ["alpha"], False, "Symmetric packed rank-2 matrix update"),
+  Routine(True,  "2b", "syr2",  T,  [S,D],     ["n"], ["layout","triangle"], ["x","y"], ["a"], ["alpha"], False, "Symmetric rank-2 matrix update"),
+  Routine(True,  "2b", "spr2",  T,  [S,D],     ["n"], ["layout","triangle"], ["x","y"], ["ap"], ["alpha"], False, "Symmetric packed rank-2 matrix update"),
 ],
 [ # Level 3: matrix-matrix
   Routine(True,  "3", "gemm",  T,  [S,D,C,Z], ["m","n","k"], ["layout","a_transpose","b_transpose"], ["a","b"], ["c"], ["alpha","beta"], False, "General matrix-matrix multiplication"),
@@ -247,7 +247,7 @@ files = [
   path_clblast+"/src/clblast_c.cc",
   path_clblast+"/test/wrapper_clblas.h",
 ]
-header_lines = [84, 59, 80, 24, 22]
+header_lines = [84, 63, 80, 24, 22]
 footer_lines = [6, 3, 5, 2, 6]
 
 # Checks whether the command-line arguments are valid; exists otherwise
diff --git a/src/clblast.cc b/src/clblast.cc
index 466de83e..3695aa02 100644
--- a/src/clblast.cc
+++ b/src/clblast.cc
@@ -43,8 +43,12 @@
 #include "internal/routines/level2/xgerc.h"
 #include "internal/routines/level2/xher.h"
 #include "internal/routines/level2/xhpr.h"
+#include "internal/routines/level2/xher2.h"
+#include "internal/routines/level2/xhpr2.h"
 #include "internal/routines/level2/xsyr.h"
 #include "internal/routines/level2/xspr.h"
+#include "internal/routines/level2/xsyr2.h"
+#include "internal/routines/level2/xspr2.h"
 
 // BLAS level-3 includes
 #include "internal/routines/level3/xgemm.h"
@@ -1014,14 +1018,24 @@ template StatusCode Hpr<double>(const Layout, const Triangle,
 
 // Hermitian rank-2 matrix update: CHER2/ZHER2
 template <typename T>
-StatusCode Her2(const Layout, const Triangle,
-                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 Her2(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,
+                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 = Xher2<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoHer2(layout, triangle,
+                        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 Her2<float2>(const Layout, const Triangle,
                                  const size_t,
@@ -1040,14 +1054,24 @@ template StatusCode Her2<double2>(const Layout, const Triangle,
 
 // Hermitian packed rank-2 matrix update: CHPR2/ZHPR2
 template <typename T>
-StatusCode Hpr2(const Layout, const Triangle,
-                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,
-                cl_command_queue*, cl_event*) {
-  return StatusCode::kNotImplemented;
+StatusCode Hpr2(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,
+                const cl_mem y_buffer, const size_t y_offset, const size_t y_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 = Xhpr2<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoHpr2(layout, triangle,
+                        n,
+                        alpha,
+                        Buffer<T>(x_buffer), x_offset, x_inc,
+                        Buffer<T>(y_buffer), y_offset, y_inc,
+                        Buffer<T>(ap_buffer), ap_offset);
 }
 template StatusCode Hpr2<float2>(const Layout, const Triangle,
                                  const size_t,
@@ -1130,14 +1154,24 @@ template StatusCode Spr<double>(const Layout, const Triangle,
 
 // Symmetric rank-2 matrix update: SSYR2/DSYR2
 template <typename T>
-StatusCode Syr2(const Layout, const Triangle,
-                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 Syr2(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,
+                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 = Xsyr2<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoSyr2(layout, triangle,
+                        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 Syr2<float>(const Layout, const Triangle,
                                 const size_t,
@@ -1156,14 +1190,24 @@ template StatusCode Syr2<double>(const Layout, const Triangle,
 
 // Symmetric packed rank-2 matrix update: SSPR2/DSPR2
 template <typename T>
-StatusCode Spr2(const Layout, const Triangle,
-                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,
-                cl_command_queue*, cl_event*) {
-  return StatusCode::kNotImplemented;
+StatusCode Spr2(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,
+                const cl_mem y_buffer, const size_t y_offset, const size_t y_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 = Xspr2<T>(queue_cpp, event_cpp);
+  auto status = routine.SetUp();
+  if (status != StatusCode::kSuccess) { return status; }
+  return routine.DoSpr2(layout, triangle,
+                        n,
+                        alpha,
+                        Buffer<T>(x_buffer), x_offset, x_inc,
+                        Buffer<T>(y_buffer), y_offset, y_inc,
+                        Buffer<T>(ap_buffer), ap_offset);
 }
 template StatusCode Spr2<float>(const Layout, const Triangle,
                                 const size_t,
diff --git a/src/kernels/common.opencl b/src/kernels/common.opencl
index 973c123e..f2a2e7a7 100644
--- a/src/kernels/common.opencl
+++ b/src/kernels/common.opencl
@@ -147,13 +147,6 @@ 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/level2.opencl b/src/kernels/level2/level2.opencl
index ad92595a..be979766 100644
--- a/src/kernels/level2/level2.opencl
+++ b/src/kernels/level2/level2.opencl
@@ -39,10 +39,7 @@ inline real LoadVector(const int id, const int max,
   if (id < max) {
     real result = gm[id*inc + offset];
     if (do_conjugate) {
-      #if defined(ROUTINE_GERC)
-        COMPLEX_CONJUGATE(result);
-      #endif
-      #if defined(ROUTINE_HER) || defined(ROUTINE_HPR)
+      #if defined(ROUTINE_GERC) || defined(ROUTINE_HER) || defined(ROUTINE_HPR) || defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
         COMPLEX_CONJUGATE(result);
       #endif
     }
@@ -81,8 +78,16 @@ inline void MatrixUpdate(const int id1, const int id2, const int max1, const int
     const real avalue = agm[a_index];
 
     // Computes result = alpha * x[i] * y[j] + a[i][j]
-    real result;
-    GER(result, alpha, xvalue, yvalue, avalue);
+    #if PRECISION == 3232 || PRECISION == 6464
+      real ax;
+      ax.x = MulReal(alpha, xvalue);
+      ax.y = MulImag(alpha, xvalue);
+      real result;
+      result.x = MulReal(ax, yvalue) + avalue.x;
+      result.y = MulImag(ax, yvalue) + avalue.y;
+    #else
+      real result = alpha * xvalue * yvalue + avalue;
+    #endif
 
     // For hermetian matrices
     #if defined(ROUTINE_HER) || defined(ROUTINE_HPR)
@@ -94,6 +99,57 @@ inline void MatrixUpdate(const int id1, const int id2, const int max1, const int
   }
 }
 
+// Performs the rank-2 matrix update
+inline void MatrixUpdate2(const int id1, const int id2, const int max1, const int max2,
+                          __global real* agm, const int a_offset, const int a_ld,
+                          const real alpha1, const real xvalue, const real yvalue,
+                          const real alpha2, const real xtvalue, const real ytvalue,
+                          const int is_upper) {
+
+  // Bounds of a regular matrix
+  if (id1 < max1 && id2 < max2) {
+
+    #if defined(ROUTINE_SPR2) || defined(ROUTINE_HPR2)
+      int a_index;
+      if (is_upper) {
+        a_index = (id1 <= id2) ? ((id2+1)*id2)/2 + id1 : ((id1+1)*id1)/2 + id2;
+      }
+      else {
+        a_index = (id1 >= id2) ? ((2*a_ld-(id2+1))*id2)/2 + id1 : ((2*a_ld-(id1+1))*id1)/2 + id2;
+      }
+      a_index += 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] + alpha * x[j] * y[i] + a[i][j]
+    #if PRECISION == 3232 || PRECISION == 6464
+      real ax;
+      ax.x = MulReal(alpha2, xvalue);
+      ax.y = MulImag(alpha2, xvalue);
+      real atx;
+      atx.x = MulReal(alpha1, xtvalue);
+      atx.y = MulImag(alpha1, xtvalue);
+      real result;
+      result.x = MulReal(ax, yvalue) + MulReal(atx, ytvalue) + avalue.x;
+      result.y = MulImag(ax, yvalue) + MulImag(atx, ytvalue) + avalue.y;
+    #else
+      real result = alpha1 * xvalue * yvalue + alpha2 * xtvalue * ytvalue + avalue;
+    #endif
+
+    // For hermetian matrices
+    #if defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
+      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/xher2.opencl b/src/kernels/level2/xher2.opencl
new file mode 100644
index 00000000..4a2edce8
--- /dev/null
+++ b/src/kernels/level2/xher2.opencl
@@ -0,0 +1,104 @@
+
+// =================================================================================================
+// 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 Xher2 kernels for rank-2 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-2 matrix update kernel (HER2, HPR2, SYR2, SPR2)
+__attribute__((reqd_work_group_size(WGS1, WGS2, 1)))
+__kernel void Xher2(const int n, const real alpha,
+                    const __global real* restrict xgm, const int x_offset, const int x_inc,
+                    const __global real* restrict ygm, const int y_offset, const int y_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 Y
+  real xvalues[WPT];
+  real yvalues[WPT];
+  real xtvalues[WPT];
+  real ytvalues[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);
+  }
+
+  // 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);
+    yvalues[w] = LoadVector(id1, n, ygm, y_offset, y_inc, is_rowmajor);
+  }
+
+  // Loads the Y-transposed-vector
+  #pragma unroll
+  for (int w=0; w<WPT; ++w) {
+    const int id2 = w*get_global_size(1) + get_global_id(1);
+    ytvalues[w] = LoadVector(id2, n, ygm, y_offset, y_inc, !is_rowmajor);
+  }
+
+  // Sets the proper value of alpha in case conjugation is needed
+  real alpha1 = alpha;
+  real alpha2 = alpha;
+  #if defined(ROUTINE_HER2) || defined(ROUTINE_HPR2)
+    if (is_rowmajor) {
+      COMPLEX_CONJUGATE(alpha1);
+    }
+    else {
+      COMPLEX_CONJUGATE(alpha2);
+    }
+  #endif
+
+  // 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 {
+        MatrixUpdate2(id1, id2, n, n, agm, a_offset, a_ld,
+                      alpha1, xvalues[w2], yvalues[w1],
+                      alpha2, xtvalues[w1], ytvalues[w2], is_upper);
+      }
+    }
+  }
+}
+
+// =================================================================================================
+
+// End of the C++11 raw string literal
+)"
+
+// =================================================================================================
diff --git a/src/routines/level2/xher2.cc b/src/routines/level2/xher2.cc
new file mode 100644
index 00000000..63144f77
--- /dev/null
+++ b/src/routines/level2/xher2.cc
@@ -0,0 +1,114 @@
+
+// =================================================================================================
+// 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 Xher2 class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xher2.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Specific implementations to get the memory-type based on a template argument
+template <> const Precision Xher2<float>::precision_ = Precision::kSingle;
+template <> const Precision Xher2<double>::precision_ = Precision::kDouble;
+template <> const Precision Xher2<float2>::precision_ = Precision::kComplexSingle;
+template <> const Precision Xher2<double2>::precision_ = Precision::kComplexDouble;
+
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xher2<T>::Xher2(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/xher2.opencl"
+  ;
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xher2<T>::DoHer2(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> &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,
+                            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);
+
+  // 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; }
+  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, "Xher2");
+
+    // Sets the kernel arguments
+    kernel.SetArgument(0, static_cast<int>(n));
+    kernel.SetArgument(1, 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, y_buffer());
+    kernel.SetArgument(6, static_cast<int>(y_offset));
+    kernel.SetArgument(7, static_cast<int>(y_inc));
+    kernel.SetArgument(8, a_buffer());
+    kernel.SetArgument(9, static_cast<int>(a_offset));
+    kernel.SetArgument(10, static_cast<int>(a_ld));
+    kernel.SetArgument(11, static_cast<int>(is_upper));
+    kernel.SetArgument(12, 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 Xher2<float>;
+template class Xher2<double>;
+template class Xher2<float2>;
+template class Xher2<double2>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xhpr.cc b/src/routines/level2/xhpr.cc
index abe00669..b0cea72f 100644
--- a/src/routines/level2/xhpr.cc
+++ b/src/routines/level2/xhpr.cc
@@ -34,10 +34,11 @@ StatusCode Xhpr<T,U>::DoHpr(const Layout layout, const Triangle triangle,
                             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) {
 
-  // 
+  // Specific Xhpr functionality is implemented in the kernel using defines
   return DoHer(layout, triangle, n, alpha,
                x_buffer, x_offset, x_inc,
-               ap_buffer, ap_offset, n, true);
+               ap_buffer, ap_offset, n,
+               true); // packed matrix
 }
 
 // =================================================================================================
diff --git a/src/routines/level2/xhpr2.cc b/src/routines/level2/xhpr2.cc
new file mode 100644
index 00000000..ded35e53
--- /dev/null
+++ b/src/routines/level2/xhpr2.cc
@@ -0,0 +1,53 @@
+
+// =================================================================================================
+// 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 Xhpr2 class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xhpr2.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xhpr2<T>::Xhpr2(Queue &queue, Event &event, const std::string &name):
+    Xher2<T>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xhpr2<T>::DoHpr2(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> &y_buffer, const size_t y_offset, const size_t y_inc,
+                            const Buffer<T> &ap_buffer, const size_t ap_offset) {
+
+  // Specific Xhpr2 functionality is implemented in the kernel using defines
+  return DoHer2(layout, triangle, n, alpha,
+                x_buffer, x_offset, x_inc,
+                y_buffer, y_offset, y_inc,
+                ap_buffer, ap_offset, n,
+                true); // packed matrix
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xhpr2<float2>;
+template class Xhpr2<double2>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xspr.cc b/src/routines/level2/xspr.cc
index 5159ad50..2d998e0b 100644
--- a/src/routines/level2/xspr.cc
+++ b/src/routines/level2/xspr.cc
@@ -34,10 +34,11 @@ StatusCode Xspr<T>::DoSpr(const Layout layout, const Triangle triangle,
                           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) {
 
-  // 
+  // Specific Xspr functionality is implemented in the kernel using defines
   return DoHer(layout, triangle, n, alpha,
                x_buffer, x_offset, x_inc,
-               ap_buffer, ap_offset, n, true);
+               ap_buffer, ap_offset, n,
+               true); // packed matrix
 }
 
 // =================================================================================================
diff --git a/src/routines/level2/xspr2.cc b/src/routines/level2/xspr2.cc
new file mode 100644
index 00000000..fd5232da
--- /dev/null
+++ b/src/routines/level2/xspr2.cc
@@ -0,0 +1,53 @@
+
+// =================================================================================================
+// 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 Xspr2 class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xspr2.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xspr2<T>::Xspr2(Queue &queue, Event &event, const std::string &name):
+    Xher2<T>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xspr2<T>::DoSpr2(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> &y_buffer, const size_t y_offset, const size_t y_inc,
+                            const Buffer<T> &ap_buffer, const size_t ap_offset) {
+
+  // Specific Xspr2 functionality is implemented in the kernel using defines
+  return DoHer2(layout, triangle, n, alpha,
+                x_buffer, x_offset, x_inc,
+                y_buffer, y_offset, y_inc,
+                ap_buffer, ap_offset, n,
+                true); // packed matrix
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xspr2<float>;
+template class Xspr2<double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/src/routines/level2/xsyr.cc b/src/routines/level2/xsyr.cc
index 755fde0d..c01fa2d3 100644
--- a/src/routines/level2/xsyr.cc
+++ b/src/routines/level2/xsyr.cc
@@ -34,7 +34,7 @@ StatusCode Xsyr<T>::DoSyr(const Layout layout, const Triangle triangle,
                           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) {
 
-  // 
+  // Specific Xsyr functionality is implemented in the kernel using defines
   return DoHer(layout, triangle, n, alpha,
                x_buffer, x_offset, x_inc,
                a_buffer, a_offset, a_ld);
diff --git a/src/routines/level2/xsyr2.cc b/src/routines/level2/xsyr2.cc
new file mode 100644
index 00000000..6db55085
--- /dev/null
+++ b/src/routines/level2/xsyr2.cc
@@ -0,0 +1,52 @@
+
+// =================================================================================================
+// 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 Xsyr2 class (see the header for information about the class).
+//
+// =================================================================================================
+
+#include "internal/routines/level2/xsyr2.h"
+
+#include <string>
+
+namespace clblast {
+// =================================================================================================
+
+// Constructor: forwards to base class constructor
+template <typename T>
+Xsyr2<T>::Xsyr2(Queue &queue, Event &event, const std::string &name):
+    Xher2<T>(queue, event, name) {
+}
+
+// =================================================================================================
+
+// The main routine
+template <typename T>
+StatusCode Xsyr2<T>::DoSyr2(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> &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) {
+
+  // Specific Xsyr2 functionality is implemented in the kernel using defines
+  return DoHer2(layout, triangle, n, alpha,
+                x_buffer, x_offset, x_inc,
+                y_buffer, y_offset, y_inc,
+                a_buffer, a_offset, a_ld);
+}
+
+// =================================================================================================
+
+// Compiles the templated class
+template class Xsyr2<float>;
+template class Xsyr2<double>;
+
+// =================================================================================================
+} // namespace clblast
diff --git a/test/routines/level2/xher2.h b/test/routines/level2/xher2.h
new file mode 100644
index 00000000..c12ff827
--- /dev/null
+++ b/test/routines/level2/xher2.h
@@ -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 implements a class with static methods to describe the Xher2 routine. Examples of
+// such 'descriptions' are how to calculate the size a of buffer or how to run the routine. These
+// static methods are used by the correctness tester and the performance tester.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_TEST_ROUTINES_XHER2_H_
+#define CLBLAST_TEST_ROUTINES_XHER2_H_
+
+#include <vector>
+#include <string>
+
+#include "wrapper_clblas.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class TestXher2 {
+ public:
+
+  // The BLAS level: 1, 2, or 3
+  static size_t BLASLevel() { return 2; }
+
+  // The list of arguments relevant for this routine
+  static std::vector<std::string> GetOptions() {
+    return {kArgN,
+            kArgLayout, kArgTriangle,
+            kArgALeadDim, kArgXInc, kArgYInc,
+            kArgAOffset, kArgXOffset, kArgYOffset,
+            kArgAlpha};
+  }
+
+  // Describes how to obtain the sizes of the buffers
+  static size_t GetSizeX(const Arguments<T> &args) {
+    return args.n * args.x_inc + args.x_offset;
+  }
+  static size_t GetSizeY(const Arguments<T> &args) {
+    return args.n * args.y_inc + args.y_offset;
+  }
+  static size_t GetSizeA(const Arguments<T> &args) {
+    return args.n * args.a_ld + args.a_offset;
+  }
+
+  // Describes how to set the sizes of all the buffers
+  static void SetSizes(Arguments<T> &args) {
+    args.a_size = GetSizeA(args);
+    args.x_size = GetSizeX(args);
+    args.y_size = GetSizeY(args);
+  }
+
+  // Describes what the default values of the leading dimensions of the matrices are
+  static size_t DefaultLDA(const Arguments<T> &args) { return args.n; }
+  static size_t DefaultLDB(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t DefaultLDC(const Arguments<T> &) { return 1; } // N/A for this routine
+
+  // Describes which transpose options are relevant for this routine
+  using Transposes = std::vector<Transpose>;
+  static Transposes GetATransposes(const Transposes &) { return {}; } // N/A for this routine
+  static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
+
+  // Describes how to run the CLBlast routine
+  static StatusCode RunRoutine(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = Her2(args.layout, args.triangle,
+                       args.n, args.alpha,
+                       buffers.x_vec(), args.x_offset, args.x_inc,
+                       buffers.y_vec(), args.y_offset, args.y_inc,
+                       buffers.a_mat(), args.a_offset, args.a_ld,
+                       &queue_plain, &event);
+    clWaitForEvents(1, &event);
+    return status;
+  }
+
+  // Describes how to run the clBLAS routine (for correctness/performance comparison)
+  static StatusCode RunReference(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXher2(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasUplo>(args.triangle),
+                              args.n, args.alpha,
+                              buffers.x_vec(), args.x_offset, args.x_inc,
+                              buffers.y_vec(), args.y_offset, args.y_inc,
+                              buffers.a_mat(), args.a_offset, args.a_ld,
+                              1, &queue_plain, 0, nullptr, &event);
+    clWaitForEvents(1, &event);
+    return static_cast<StatusCode>(status);
+  }
+
+  // Describes how to download the results of the computation (more importantly: which buffer)
+  static std::vector<T> DownloadResult(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+    std::vector<T> result(args.a_size, static_cast<T>(0));
+    buffers.a_mat.Read(queue, args.a_size, result);
+    return result;
+  }
+
+  // Describes how to compute the indices of the result buffer
+  static size_t ResultID1(const Arguments<T> &args) { return args.n; }
+  static size_t ResultID2(const Arguments<T> &args) { return args.n; }
+  static size_t GetResultIndex(const Arguments<T> &args, const size_t id1, const size_t id2) {
+    return id2*args.a_ld + id1 + args.a_offset;
+  }
+
+  // Describes how to compute performance metrics
+  static size_t GetFlops(const Arguments<T> &args) {
+    return 5 * args.n * args.n;
+  }
+  static size_t GetBytes(const Arguments<T> &args) {
+    return (args.n*args.n + 2 * args.n) * sizeof(T);
+  }
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_TEST_ROUTINES_XHER2_H_
+#endif
diff --git a/test/routines/level2/xhpr2.h b/test/routines/level2/xhpr2.h
new file mode 100644
index 00000000..68fbc76d
--- /dev/null
+++ b/test/routines/level2/xhpr2.h
@@ -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 implements a class with static methods to describe the Xhpr2 routine. Examples of
+// such 'descriptions' are how to calculate the size a of buffer or how to run the routine. These
+// static methods are used by the correctness tester and the performance tester.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_TEST_ROUTINES_XHPR2_H_
+#define CLBLAST_TEST_ROUTINES_XHPR2_H_
+
+#include <vector>
+#include <string>
+
+#include "wrapper_clblas.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class TestXhpr2 {
+ public:
+
+  // The BLAS level: 1, 2, or 3
+  static size_t BLASLevel() { return 2; }
+
+  // The list of arguments relevant for this routine
+  static std::vector<std::string> GetOptions() {
+    return {kArgN,
+            kArgLayout, kArgTriangle,
+            kArgXInc, kArgYInc,
+            kArgAPOffset, kArgXOffset, kArgYOffset,
+            kArgAlpha};
+  }
+
+  // Describes how to obtain the sizes of the buffers
+  static size_t GetSizeX(const Arguments<T> &args) {
+    return args.n * args.x_inc + args.x_offset;
+  }
+  static size_t GetSizeY(const Arguments<T> &args) {
+    return args.n * args.y_inc + args.y_offset;
+  }
+  static size_t GetSizeAP(const Arguments<T> &args) {
+    return ((args.n*(args.n+1)) / 2) + args.ap_offset;
+  }
+
+  // Describes how to set the sizes of all the buffers
+  static void SetSizes(Arguments<T> &args) {
+    args.ap_size = GetSizeAP(args);
+    args.x_size = GetSizeX(args);
+    args.y_size = GetSizeY(args);
+  }
+
+  // Describes what the default values of the leading dimensions of the matrices are
+  static size_t DefaultLDA(const Arguments<T> &args) { return args.n; }
+  static size_t DefaultLDB(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t DefaultLDC(const Arguments<T> &) { return 1; } // N/A for this routine
+
+  // Describes which transpose options are relevant for this routine
+  using Transposes = std::vector<Transpose>;
+  static Transposes GetATransposes(const Transposes &) { return {}; } // N/A for this routine
+  static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
+
+  // Describes how to run the CLBlast routine
+  static StatusCode RunRoutine(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = Hpr2(args.layout, args.triangle,
+                       args.n, args.alpha,
+                       buffers.x_vec(), args.x_offset, args.x_inc,
+                       buffers.y_vec(), args.y_offset, args.y_inc,
+                       buffers.ap_mat(), args.ap_offset,
+                       &queue_plain, &event);
+    clWaitForEvents(1, &event);
+    return status;
+  }
+
+  // Describes how to run the clBLAS routine (for correctness/performance comparison)
+  static StatusCode RunReference(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXhpr2(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasUplo>(args.triangle),
+                              args.n, args.alpha,
+                              buffers.x_vec(), args.x_offset, args.x_inc,
+                              buffers.y_vec(), args.y_offset, args.y_inc,
+                              buffers.ap_mat(), args.ap_offset,
+                              1, &queue_plain, 0, nullptr, &event);
+    clWaitForEvents(1, &event);
+    return static_cast<StatusCode>(status);
+  }
+
+  // Describes how to download the results of the computation (more importantly: which buffer)
+  static std::vector<T> DownloadResult(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+    std::vector<T> result(args.ap_size, static_cast<T>(0));
+    buffers.ap_mat.Read(queue, args.ap_size, result);
+    return result;
+  }
+
+  // Describes how to compute the indices of the result buffer
+  static size_t ResultID1(const Arguments<T> &args) { return args.ap_size - args.ap_offset; }
+  static size_t ResultID2(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t GetResultIndex(const Arguments<T> &args, const size_t id1, const size_t) {
+    return id1 + args.ap_offset;
+  }
+
+  // Describes how to compute performance metrics
+  static size_t GetFlops(const Arguments<T> &args) {
+    return 5 * ((args.n*(args.n+1)) / 2);
+  }
+  static size_t GetBytes(const Arguments<T> &args) {
+    return ((args.n*(args.n+1)) + 2 * args.n) * sizeof(T);
+  }
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_TEST_ROUTINES_XHPR2_H_
+#endif
diff --git a/test/routines/level2/xspr2.h b/test/routines/level2/xspr2.h
new file mode 100644
index 00000000..43d66c9e
--- /dev/null
+++ b/test/routines/level2/xspr2.h
@@ -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 implements a class with static methods to describe the Xspr2 routine. Examples of
+// such 'descriptions' are how to calculate the size a of buffer or how to run the routine. These
+// static methods are used by the correctness tester and the performance tester.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_TEST_ROUTINES_XSPR2_H_
+#define CLBLAST_TEST_ROUTINES_XSPR2_H_
+
+#include <vector>
+#include <string>
+
+#include "wrapper_clblas.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class TestXspr2 {
+ public:
+
+  // The BLAS level: 1, 2, or 3
+  static size_t BLASLevel() { return 2; }
+
+  // The list of arguments relevant for this routine
+  static std::vector<std::string> GetOptions() {
+    return {kArgN,
+            kArgLayout, kArgTriangle,
+            kArgXInc, kArgYInc,
+            kArgAPOffset, kArgXOffset, kArgYOffset,
+            kArgAlpha};
+  }
+
+  // Describes how to obtain the sizes of the buffers
+  static size_t GetSizeX(const Arguments<T> &args) {
+    return args.n * args.x_inc + args.x_offset;
+  }
+  static size_t GetSizeY(const Arguments<T> &args) {
+    return args.n * args.y_inc + args.y_offset;
+  }
+  static size_t GetSizeAP(const Arguments<T> &args) {
+    return ((args.n*(args.n+1)) / 2) + args.ap_offset;
+  }
+
+  // Describes how to set the sizes of all the buffers
+  static void SetSizes(Arguments<T> &args) {
+    args.ap_size = GetSizeAP(args);
+    args.x_size = GetSizeX(args);
+    args.y_size = GetSizeY(args);
+  }
+
+  // Describes what the default values of the leading dimensions of the matrices are
+  static size_t DefaultLDA(const Arguments<T> &args) { return args.n; }
+  static size_t DefaultLDB(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t DefaultLDC(const Arguments<T> &) { return 1; } // N/A for this routine
+
+  // Describes which transpose options are relevant for this routine
+  using Transposes = std::vector<Transpose>;
+  static Transposes GetATransposes(const Transposes &) { return {}; } // N/A for this routine
+  static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
+
+  // Describes how to run the CLBlast routine
+  static StatusCode RunRoutine(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = Spr2(args.layout, args.triangle,
+                       args.n, args.alpha,
+                       buffers.x_vec(), args.x_offset, args.x_inc,
+                       buffers.y_vec(), args.y_offset, args.y_inc,
+                       buffers.ap_mat(), args.ap_offset,
+                       &queue_plain, &event);
+    clWaitForEvents(1, &event);
+    return status;
+  }
+
+  // Describes how to run the clBLAS routine (for correctness/performance comparison)
+  static StatusCode RunReference(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXspr2(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasUplo>(args.triangle),
+                              args.n, args.alpha,
+                              buffers.x_vec(), args.x_offset, args.x_inc,
+                              buffers.y_vec(), args.y_offset, args.y_inc,
+                              buffers.ap_mat(), args.ap_offset,
+                              1, &queue_plain, 0, nullptr, &event);
+    clWaitForEvents(1, &event);
+    return static_cast<StatusCode>(status);
+  }
+
+  // Describes how to download the results of the computation (more importantly: which buffer)
+  static std::vector<T> DownloadResult(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+    std::vector<T> result(args.ap_size, static_cast<T>(0));
+    buffers.ap_mat.Read(queue, args.ap_size, result);
+    return result;
+  }
+
+  // Describes how to compute the indices of the result buffer
+  static size_t ResultID1(const Arguments<T> &args) { return args.ap_size - args.ap_offset; }
+  static size_t ResultID2(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t GetResultIndex(const Arguments<T> &args, const size_t id1, const size_t) {
+    return id1 + args.ap_offset;
+  }
+
+  // Describes how to compute performance metrics
+  static size_t GetFlops(const Arguments<T> &args) {
+    return 5 * ((args.n*(args.n+1)) / 2);
+  }
+  static size_t GetBytes(const Arguments<T> &args) {
+    return ((args.n*(args.n+1)) + 2 * args.n) * sizeof(T);
+  }
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_TEST_ROUTINES_XSPR2_H_
+#endif
diff --git a/test/routines/level2/xsyr2.h b/test/routines/level2/xsyr2.h
new file mode 100644
index 00000000..32497a61
--- /dev/null
+++ b/test/routines/level2/xsyr2.h
@@ -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 implements a class with static methods to describe the Xsyr2 routine. Examples of
+// such 'descriptions' are how to calculate the size a of buffer or how to run the routine. These
+// static methods are used by the correctness tester and the performance tester.
+//
+// =================================================================================================
+
+#ifndef CLBLAST_TEST_ROUTINES_XSYR2_H_
+#define CLBLAST_TEST_ROUTINES_XSYR2_H_
+
+#include <vector>
+#include <string>
+
+#include "wrapper_clblas.h"
+
+namespace clblast {
+// =================================================================================================
+
+// See comment at top of file for a description of the class
+template <typename T>
+class TestXsyr2 {
+ public:
+
+  // The BLAS level: 1, 2, or 3
+  static size_t BLASLevel() { return 2; }
+
+  // The list of arguments relevant for this routine
+  static std::vector<std::string> GetOptions() {
+    return {kArgN,
+            kArgLayout, kArgTriangle,
+            kArgALeadDim, kArgXInc, kArgYInc,
+            kArgAOffset, kArgXOffset, kArgYOffset,
+            kArgAlpha};
+  }
+
+  // Describes how to obtain the sizes of the buffers
+  static size_t GetSizeX(const Arguments<T> &args) {
+    return args.n * args.x_inc + args.x_offset;
+  }
+  static size_t GetSizeY(const Arguments<T> &args) {
+    return args.n * args.y_inc + args.y_offset;
+  }
+  static size_t GetSizeA(const Arguments<T> &args) {
+    return args.n * args.a_ld + args.a_offset;
+  }
+
+  // Describes how to set the sizes of all the buffers
+  static void SetSizes(Arguments<T> &args) {
+    args.a_size = GetSizeA(args);
+    args.x_size = GetSizeX(args);
+    args.y_size = GetSizeY(args);
+  }
+
+  // Describes what the default values of the leading dimensions of the matrices are
+  static size_t DefaultLDA(const Arguments<T> &args) { return args.n; }
+  static size_t DefaultLDB(const Arguments<T> &) { return 1; } // N/A for this routine
+  static size_t DefaultLDC(const Arguments<T> &) { return 1; } // N/A for this routine
+
+  // Describes which transpose options are relevant for this routine
+  using Transposes = std::vector<Transpose>;
+  static Transposes GetATransposes(const Transposes &) { return {}; } // N/A for this routine
+  static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
+
+  // Describes how to run the CLBlast routine
+  static StatusCode RunRoutine(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = Syr2(args.layout, args.triangle,
+                       args.n, args.alpha,
+                       buffers.x_vec(), args.x_offset, args.x_inc,
+                       buffers.y_vec(), args.y_offset, args.y_inc,
+                       buffers.a_mat(), args.a_offset, args.a_ld,
+                       &queue_plain, &event);
+    clWaitForEvents(1, &event);
+    return status;
+  }
+
+  // Describes how to run the clBLAS routine (for correctness/performance comparison)
+  static StatusCode RunReference(const Arguments<T> &args, const Buffers<T> &buffers, Queue &queue) {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXsyr2(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasUplo>(args.triangle),
+                              args.n, args.alpha,
+                              buffers.x_vec(), args.x_offset, args.x_inc,
+                              buffers.y_vec(), args.y_offset, args.y_inc,
+                              buffers.a_mat(), args.a_offset, args.a_ld,
+                              1, &queue_plain, 0, nullptr, &event);
+    clWaitForEvents(1, &event);
+    return static_cast<StatusCode>(status);
+  }
+
+  // Describes how to download the results of the computation (more importantly: which buffer)
+  static std::vector<T> DownloadResult(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+    std::vector<T> result(args.a_size, static_cast<T>(0));
+    buffers.a_mat.Read(queue, args.a_size, result);
+    return result;
+  }
+
+  // Describes how to compute the indices of the result buffer
+  static size_t ResultID1(const Arguments<T> &args) { return args.n; }
+  static size_t ResultID2(const Arguments<T> &args) { return args.n; }
+  static size_t GetResultIndex(const Arguments<T> &args, const size_t id1, const size_t id2) {
+    return id2*args.a_ld + id1 + args.a_offset;
+  }
+
+  // Describes how to compute performance metrics
+  static size_t GetFlops(const Arguments<T> &args) {
+    return 5 * args.n * args.n;
+  }
+  static size_t GetBytes(const Arguments<T> &args) {
+    return (args.n*args.n + 2 * args.n) * sizeof(T);
+  }
+};
+
+// =================================================================================================
+} // namespace clblast
+
+// CLBLAST_TEST_ROUTINES_XSYR2_H_
+#endif