Added support for testing (performance and correctness) against a CPU BLAS library

22 files changed, 902 insertions, 326 deletions

diff --git a/test/routines/level2/xgbmv.h b/test/routines/level2/xgbmv.h
index 0e238804..b875075d 100644
--- a/test/routines/level2/xgbmv.h
+++ b/test/routines/level2/xgbmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -76,7 +81,7 @@ class TestXgbmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Gbmv(args.layout, args.a_transpose,
@@ -90,19 +95,41 @@ class TestXgbmv {
   }
 
   // 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 = clblasXgbmv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasTranspose>(args.a_transpose),
-                              args.m, args.n, args.kl, args.ku, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXgbmv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasTranspose>(args.a_transpose),
+                                args.m, args.n, args.kl, args.ku, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXgbmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.a_transpose),
+                 args.m, args.n, args.kl, args.ku, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xgemv.h b/test/routines/level2/xgemv.h
index 2924d498..a70ccd34 100644
--- a/test/routines/level2/xgemv.h
+++ b/test/routines/level2/xgemv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -76,7 +81,7 @@ class TestXgemv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Gemv(args.layout, args.a_transpose,
@@ -90,19 +95,41 @@ class TestXgemv {
   }
 
   // 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 = clblasXgemv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasTranspose>(args.a_transpose),
-                              args.m, args.n, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXgemv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasTranspose>(args.a_transpose),
+                                args.m, args.n, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXgemv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.a_transpose),
+                 args.m, args.n, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xger.h b/test/routines/level2/xger.h
index 98296e92..32c2a505 100644
--- a/test/routines/level2/xger.h
+++ b/test/routines/level2/xger.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -72,7 +77,7 @@ class TestXger {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Ger(args.layout,
@@ -86,18 +91,39 @@ class TestXger {
   }
 
   // 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 = clblasXger(static_cast<clblasOrder>(args.layout),
-                             args.m, 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXger(static_cast<clblasOrder>(args.layout),
+                               args.m, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXger(convertToCBLAS(args.layout),
+                args.m, args.n, args.alpha,
+                x_vec_cpu, args.x_offset, args.x_inc,
+                y_vec_cpu, args.y_offset, args.y_inc,
+                a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xgerc.h b/test/routines/level2/xgerc.h
index 77258d92..4b6954f6 100644
--- a/test/routines/level2/xgerc.h
+++ b/test/routines/level2/xgerc.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -72,7 +77,7 @@ class TestXgerc {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Gerc(args.layout,
@@ -86,18 +91,39 @@ class TestXgerc {
   }
 
   // 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 = clblasXgerc(static_cast<clblasOrder>(args.layout),
-                              args.m, 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXgerc(static_cast<clblasOrder>(args.layout),
+                                args.m, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXgerc(convertToCBLAS(args.layout),
+                 args.m, args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xgeru.h b/test/routines/level2/xgeru.h
index e5f5f235..295e69e5 100644
--- a/test/routines/level2/xgeru.h
+++ b/test/routines/level2/xgeru.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -72,7 +77,7 @@ class TestXgeru {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Geru(args.layout,
@@ -86,18 +91,39 @@ class TestXgeru {
   }
 
   // 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 = clblasXgeru(static_cast<clblasOrder>(args.layout),
-                              args.m, 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXgeru(static_cast<clblasOrder>(args.layout),
+                                args.m, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXgeru(convertToCBLAS(args.layout),
+                 args.m, args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xhbmv.h b/test/routines/level2/xhbmv.h
index 34e1502f..e0bdc4da 100644
--- a/test/routines/level2/xhbmv.h
+++ b/test/routines/level2/xhbmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXhbmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Hbmv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXhbmv {
   }
 
   // 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 = clblasXhbmv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.kl, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXhbmv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.kl, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXhbmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.kl, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xhemv.h b/test/routines/level2/xhemv.h
index 80e22157..fa242961 100644
--- a/test/routines/level2/xhemv.h
+++ b/test/routines/level2/xhemv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXhemv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Hemv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXhemv {
   }
 
   // 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 = clblasXhemv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXhemv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXhemv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xher.h b/test/routines/level2/xher.h
index 53c4200f..7d0e8cc3 100644
--- a/test/routines/level2/xher.h
+++ b/test/routines/level2/xher.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -66,7 +71,7 @@ class TestXher {
   static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
 
   // Describes how to run the CLBlast routine
-  static StatusCode RunRoutine(const Arguments<U> &args, const Buffers<T> &buffers, Queue &queue) {
+  static StatusCode RunRoutine(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Her(args.layout, args.triangle,
@@ -79,18 +84,37 @@ class TestXher {
   }
 
   // Describes how to run the clBLAS routine (for correctness/performance comparison)
-  static StatusCode RunReference(const Arguments<U> &args, const Buffers<T> &buffers, Queue &queue) {
-    auto queue_plain = queue();
-    auto event = cl_event{};
-    auto status = clblasXher(static_cast<clblasOrder>(args.layout),
-                             static_cast<clblasUplo>(args.triangle),
-                             args.n, args.alpha,
-                             buffers.x_vec(), args.x_offset, args.x_inc,
-                             buffers.a_mat(), args.a_offset, args.a_ld,
-                             1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXher(static_cast<clblasOrder>(args.layout),
+                               static_cast<clblasUplo>(args.triangle),
+                               args.n, args.alpha,
+                               buffers.x_vec(), args.x_offset, args.x_inc,
+                               buffers.a_mat(), args.a_offset, args.a_ld,
+                               1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXher(convertToCBLAS(args.layout),
+                convertToCBLAS(args.triangle),
+                args.n, args.alpha,
+                x_vec_cpu, args.x_offset, args.x_inc,
+                a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // Describes how to download the results of the computation (more importantly: which buffer)
   static std::vector<T> DownloadResult(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
diff --git a/test/routines/level2/xher2.h b/test/routines/level2/xher2.h
index c12ff827..445bba74 100644
--- a/test/routines/level2/xher2.h
+++ b/test/routines/level2/xher2.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXher2 {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Her2(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXher2 {
   }
 
   // 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXher2(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xhpmv.h b/test/routines/level2/xhpmv.h
index 8fd85b62..406e564f 100644
--- a/test/routines/level2/xhpmv.h
+++ b/test/routines/level2/xhpmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXhpmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Hpmv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXhpmv {
   }
 
   // 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 = clblasXhpmv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.alpha,
-                              buffers.ap_mat(), args.ap_offset,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXhpmv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.alpha,
+                                buffers.ap_mat(), args.ap_offset,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXhpmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 ap_mat_cpu, args.ap_offset,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xhpr.h b/test/routines/level2/xhpr.h
index 03599ddc..6f56d3f3 100644
--- a/test/routines/level2/xhpr.h
+++ b/test/routines/level2/xhpr.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -66,7 +71,7 @@ class TestXhpr {
   static Transposes GetBTransposes(const Transposes &) { return {}; } // N/A for this routine
 
   // Describes how to run the CLBlast routine
-  static StatusCode RunRoutine(const Arguments<U> &args, const Buffers<T> &buffers, Queue &queue) {
+  static StatusCode RunRoutine(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Hpr(args.layout, args.triangle,
@@ -79,18 +84,37 @@ class TestXhpr {
   }
 
   // Describes how to run the clBLAS routine (for correctness/performance comparison)
-  static StatusCode RunReference(const Arguments<U> &args, const Buffers<T> &buffers, Queue &queue) {
-    auto queue_plain = queue();
-    auto event = cl_event{};
-    auto status = clblasXhpr(static_cast<clblasOrder>(args.layout),
-                             static_cast<clblasUplo>(args.triangle),
-                             args.n, args.alpha,
-                             buffers.x_vec(), args.x_offset, args.x_inc,
-                             buffers.ap_mat(), args.ap_offset,
-                             1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXhpr(static_cast<clblasOrder>(args.layout),
+                               static_cast<clblasUplo>(args.triangle),
+                               args.n, args.alpha,
+                               buffers.x_vec(), args.x_offset, args.x_inc,
+                               buffers.ap_mat(), args.ap_offset,
+                               1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXhpr(convertToCBLAS(args.layout),
+                convertToCBLAS(args.triangle),
+                args.n, args.alpha,
+                x_vec_cpu, args.x_offset, args.x_inc,
+                ap_mat_cpu, args.ap_offset);
+      buffers.ap_mat.Write(queue, args.ap_size, ap_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // Describes how to download the results of the computation (more importantly: which buffer)
   static std::vector<T> DownloadResult(const Arguments<U> &args, Buffers<T> &buffers, Queue &queue) {
diff --git a/test/routines/level2/xhpr2.h b/test/routines/level2/xhpr2.h
index 68fbc76d..43889cb9 100644
--- a/test/routines/level2/xhpr2.h
+++ b/test/routines/level2/xhpr2.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXhpr2 {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Hpr2(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXhpr2 {
   }
 
   // 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXhpr2(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 ap_mat_cpu, args.ap_offset);
+      buffers.ap_mat.Write(queue, args.ap_size, ap_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xsbmv.h b/test/routines/level2/xsbmv.h
index 5bc17e49..9a5c5140 100644
--- a/test/routines/level2/xsbmv.h
+++ b/test/routines/level2/xsbmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXsbmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Sbmv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXsbmv {
   }
 
   // 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 = clblasXsbmv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.kl, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXsbmv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.kl, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXsbmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.kl, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xspmv.h b/test/routines/level2/xspmv.h
index e335da42..913af0cd 100644
--- a/test/routines/level2/xspmv.h
+++ b/test/routines/level2/xspmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXspmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Spmv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXspmv {
   }
 
   // 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 = clblasXspmv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.alpha,
-                              buffers.ap_mat(), args.ap_offset,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXspmv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.alpha,
+                                buffers.ap_mat(), args.ap_offset,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXspmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 ap_mat_cpu, args.ap_offset,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xspr.h b/test/routines/level2/xspr.h
index 819b1ca8..bab5c541 100644
--- a/test/routines/level2/xspr.h
+++ b/test/routines/level2/xspr.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -66,7 +71,7 @@ class TestXspr {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Spr(args.layout, args.triangle,
@@ -79,18 +84,37 @@ class TestXspr {
   }
 
   // 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 = clblasXspr(static_cast<clblasOrder>(args.layout),
-                             static_cast<clblasUplo>(args.triangle),
-                             args.n, args.alpha,
-                             buffers.x_vec(), args.x_offset, args.x_inc,
-                             buffers.ap_mat(), args.ap_offset,
-                             1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXspr(static_cast<clblasOrder>(args.layout),
+                               static_cast<clblasUplo>(args.triangle),
+                               args.n, args.alpha,
+                               buffers.x_vec(), args.x_offset, args.x_inc,
+                               buffers.ap_mat(), args.ap_offset,
+                               1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXspr(convertToCBLAS(args.layout),
+                convertToCBLAS(args.triangle),
+                args.n, args.alpha,
+                x_vec_cpu, args.x_offset, args.x_inc,
+                ap_mat_cpu, args.ap_offset);
+      buffers.ap_mat.Write(queue, args.ap_size, ap_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xspr2.h b/test/routines/level2/xspr2.h
index 43d66c9e..41a04cc0 100644
--- a/test/routines/level2/xspr2.h
+++ b/test/routines/level2/xspr2.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXspr2 {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Spr2(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXspr2 {
   }
 
   // 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXspr2(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 ap_mat_cpu, args.ap_offset);
+      buffers.ap_mat.Write(queue, args.ap_size, ap_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xsymv.h b/test/routines/level2/xsymv.h
index 13473a3e..0576bc1f 100644
--- a/test/routines/level2/xsymv.h
+++ b/test/routines/level2/xsymv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXsymv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Symv(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXsymv {
   }
 
   // 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 = clblasXsymv(static_cast<clblasOrder>(args.layout),
-                              static_cast<clblasUplo>(args.triangle),
-                              args.n, args.alpha,
-                              buffers.a_mat(), args.a_offset, args.a_ld,
-                              buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
-                              buffers.y_vec(), args.y_offset, args.y_inc,
-                              1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXsymv(static_cast<clblasOrder>(args.layout),
+                                static_cast<clblasUplo>(args.triangle),
+                                args.n, args.alpha,
+                                buffers.a_mat(), args.a_offset, args.a_ld,
+                                buffers.x_vec(), args.x_offset, args.x_inc, args.beta,
+                                buffers.y_vec(), args.y_offset, args.y_inc,
+                                1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXsymv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc, args.beta,
+                 y_vec_cpu, args.y_offset, args.y_inc);
+      buffers.y_vec.Write(queue, args.y_size, y_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xsyr.h b/test/routines/level2/xsyr.h
index 66b75c0c..062eea5a 100644
--- a/test/routines/level2/xsyr.h
+++ b/test/routines/level2/xsyr.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -66,7 +71,7 @@ class TestXsyr {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Syr(args.layout, args.triangle,
@@ -79,18 +84,37 @@ class TestXsyr {
   }
 
   // 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 = clblasXsyr(static_cast<clblasOrder>(args.layout),
-                             static_cast<clblasUplo>(args.triangle),
-                             args.n, args.alpha,
-                             buffers.x_vec(), args.x_offset, args.x_inc,
-                             buffers.a_mat(), args.a_offset, args.a_ld,
-                             1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXsyr(static_cast<clblasOrder>(args.layout),
+                               static_cast<clblasUplo>(args.triangle),
+                               args.n, args.alpha,
+                               buffers.x_vec(), args.x_offset, args.x_inc,
+                               buffers.a_mat(), args.a_offset, args.a_ld,
+                               1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXsyr(convertToCBLAS(args.layout),
+                convertToCBLAS(args.triangle),
+                args.n, args.alpha,
+                x_vec_cpu, args.x_offset, args.x_inc,
+                a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xsyr2.h b/test/routines/level2/xsyr2.h
index 32497a61..50bc3cea 100644
--- a/test/routines/level2/xsyr2.h
+++ b/test/routines/level2/xsyr2.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -70,7 +75,7 @@ class TestXsyr2 {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Syr2(args.layout, args.triangle,
@@ -84,19 +89,41 @@ class TestXsyr2 {
   }
 
   // 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);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, 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);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      std::vector<T> y_vec_cpu(args.y_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      buffers.y_vec.Read(queue, args.y_size, y_vec_cpu);
+      cblasXsyr2(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 args.n, args.alpha,
+                 x_vec_cpu, args.x_offset, args.x_inc,
+                 y_vec_cpu, args.y_offset, args.y_inc,
+                 a_mat_cpu, args.a_offset, args.a_ld);
+      buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xtbmv.h b/test/routines/level2/xtbmv.h
index dbdddb65..600b4131 100644
--- a/test/routines/level2/xtbmv.h
+++ b/test/routines/level2/xtbmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -65,7 +70,7 @@ class TestXtbmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Tbmv<T>(args.layout, args.triangle, args.a_transpose, args.diagonal,
@@ -78,20 +83,41 @@ class TestXtbmv {
   }
 
   // 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 = clblasXtbmv<T>(static_cast<clblasOrder>(args.layout),
-                                 static_cast<clblasUplo>(args.triangle),
-                                 static_cast<clblasTranspose>(args.a_transpose),
-                                 static_cast<clblasDiag>(args.diagonal),
-                                 args.n, args.kl,
-                                 buffers.a_mat(), args.a_offset, args.a_ld,
-                                 buffers.x_vec(), args.x_offset, args.x_inc,
-                                 1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXtbmv<T>(static_cast<clblasOrder>(args.layout),
+                                   static_cast<clblasUplo>(args.triangle),
+                                   static_cast<clblasTranspose>(args.a_transpose),
+                                   static_cast<clblasDiag>(args.diagonal),
+                                   args.n, args.kl,
+                                   buffers.a_mat(), args.a_offset, args.a_ld,
+                                   buffers.x_vec(), args.x_offset, args.x_inc,
+                                   1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXtbmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 convertToCBLAS(args.a_transpose),
+                 convertToCBLAS(args.diagonal),
+                 args.n, args.kl,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc);
+      buffers.x_vec.Write(queue, args.x_size, x_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xtpmv.h b/test/routines/level2/xtpmv.h
index 4425765e..fc0cf393 100644
--- a/test/routines/level2/xtpmv.h
+++ b/test/routines/level2/xtpmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -65,7 +70,7 @@ class TestXtpmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Tpmv<T>(args.layout, args.triangle, args.a_transpose, args.diagonal,
@@ -78,20 +83,41 @@ class TestXtpmv {
   }
 
   // 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 = clblasXtpmv<T>(static_cast<clblasOrder>(args.layout),
-                                 static_cast<clblasUplo>(args.triangle),
-                                 static_cast<clblasTranspose>(args.a_transpose),
-                                 static_cast<clblasDiag>(args.diagonal),
-                                 args.n,
-                                 buffers.ap_mat(), args.ap_offset,
-                                 buffers.x_vec(), args.x_offset, args.x_inc,
-                                 1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXtpmv<T>(static_cast<clblasOrder>(args.layout),
+                                   static_cast<clblasUplo>(args.triangle),
+                                   static_cast<clblasTranspose>(args.a_transpose),
+                                   static_cast<clblasDiag>(args.diagonal),
+                                   args.n,
+                                   buffers.ap_mat(), args.ap_offset,
+                                   buffers.x_vec(), args.x_offset, args.x_inc,
+                                   1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> ap_mat_cpu(args.ap_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.ap_mat.Read(queue, args.ap_size, ap_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXtpmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 convertToCBLAS(args.a_transpose),
+                 convertToCBLAS(args.diagonal),
+                 args.n,
+                 ap_mat_cpu, args.ap_offset,
+                 x_vec_cpu, args.x_offset, args.x_inc);
+      buffers.x_vec.Write(queue, args.x_size, x_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {
diff --git a/test/routines/level2/xtrmv.h b/test/routines/level2/xtrmv.h
index 1c0c6fd8..fec72124 100644
--- a/test/routines/level2/xtrmv.h
+++ b/test/routines/level2/xtrmv.h
@@ -19,7 +19,12 @@
 #include <vector>
 #include <string>
 
-#include "wrapper_clblas.h"
+#ifdef CLBLAST_REF_CLBLAS
+  #include "wrapper_clblas.h"
+#endif
+#ifdef CLBLAST_REF_CBLAS
+  #include "wrapper_cblas.h"
+#endif
 
 namespace clblast {
 // =================================================================================================
@@ -65,7 +70,7 @@ class TestXtrmv {
   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) {
+  static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Trmv<T>(args.layout, args.triangle, args.a_transpose, args.diagonal,
@@ -78,20 +83,41 @@ class TestXtrmv {
   }
 
   // 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 = clblasXtrmv<T>(static_cast<clblasOrder>(args.layout),
-                                 static_cast<clblasUplo>(args.triangle),
-                                 static_cast<clblasTranspose>(args.a_transpose),
-                                 static_cast<clblasDiag>(args.diagonal),
-                                 args.n,
-                                 buffers.a_mat(), args.a_offset, args.a_ld,
-                                 buffers.x_vec(), args.x_offset, args.x_inc,
-                                 1, &queue_plain, 0, nullptr, &event);
-    clWaitForEvents(1, &event);
-    return static_cast<StatusCode>(status);
-  }
+  #ifdef CLBLAST_REF_CLBLAS
+    static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      auto queue_plain = queue();
+      auto event = cl_event{};
+      auto status = clblasXtrmv<T>(static_cast<clblasOrder>(args.layout),
+                                   static_cast<clblasUplo>(args.triangle),
+                                   static_cast<clblasTranspose>(args.a_transpose),
+                                   static_cast<clblasDiag>(args.diagonal),
+                                   args.n,
+                                   buffers.a_mat(), args.a_offset, args.a_ld,
+                                   buffers.x_vec(), args.x_offset, args.x_inc,
+                                   1, &queue_plain, 0, nullptr, &event);
+      clWaitForEvents(1, &event);
+      return static_cast<StatusCode>(status);
+    }
+  #endif
+
+  // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
+  #ifdef CLBLAST_REF_CBLAS
+    static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+      std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+      buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+      buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+      cblasXtrmv(convertToCBLAS(args.layout),
+                 convertToCBLAS(args.triangle),
+                 convertToCBLAS(args.a_transpose),
+                 convertToCBLAS(args.diagonal),
+                 args.n,
+                 a_mat_cpu, args.a_offset, args.a_ld,
+                 x_vec_cpu, args.x_offset, args.x_inc);
+      buffers.x_vec.Write(queue, args.x_size, x_vec_cpu);
+      return StatusCode::kSuccess;
+    }
+  #endif
 
   // 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) {