Merge pull request #117 from intelfx/exceptions

Convert to use C++ exceptions internally

16 files changed, 803 insertions, 923 deletions

diff --git a/src/routines/level3/xgemm.cpp b/src/routines/level3/xgemm.cpp
index 1602c69f..4f70dc7a 100644
--- a/src/routines/level3/xgemm.cpp
+++ b/src/routines/level3/xgemm.cpp
@@ -24,8 +24,7 @@ template <typename T>
 Xgemm<T>::Xgemm(Queue &queue, EventPointer event, const std::string &name):
     Routine(queue, event, name,
             {"Copy","Pad","Transpose","Padtranspose","Xgemm","XgemmDirect","KernelSelection"},
-            PrecisionValue<T>()) {
-  source_string_ =
+            PrecisionValue<T>(), {}, {
     #include "../../kernels/level3/level3.opencl"
     #include "../../kernels/level3/copy_fast.opencl"
     #include "../../kernels/level3/copy_pad.opencl"
@@ -37,30 +36,28 @@ Xgemm<T>::Xgemm(Queue &queue, EventPointer event, const std::string &name):
     #include "../../kernels/level3/xgemm_direct_part1.opencl"
     #include "../../kernels/level3/xgemm_direct_part2.opencl"
     #include "../../kernels/level3/xgemm_direct_part3.opencl"
-  ;
-  auto source_string_part_2 = // separated in two parts to prevent C1091 in MSVC 2013
+    , // separated in two parts to prevent C1091 in MSVC 2013
     #include "../../kernels/level3/xgemm_part1.opencl"
     #include "../../kernels/level3/xgemm_part2.opencl"
     #include "../../kernels/level3/xgemm_part3.opencl"
-  ;
-  source_string_ += source_string_part_2;
+    }) {
 }
 
 // =================================================================================================
 
 // The main routine
 template <typename T>
-StatusCode Xgemm<T>::DoGemm(const Layout layout,
-                            const Transpose a_transpose, const Transpose b_transpose,
-                            const size_t m, const size_t n, const size_t k,
-                            const T alpha,
-                            const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                            const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                            const T beta,
-                            const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
+void Xgemm<T>::DoGemm(const Layout layout,
+                      const Transpose a_transpose, const Transpose b_transpose,
+                      const size_t m, const size_t n, const size_t k,
+                      const T alpha,
+                      const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+                      const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+                      const T beta,
+                      const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((m == 0) || (n == 0) || (k == 0)) { return StatusCode::kInvalidDimension; }
+  if ((m == 0) || (n == 0) || (k == 0)) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes whether or not the matrices are transposed in memory. This is based on their layout
   // (row or column-major) and whether or not they are requested to be pre-transposed. Note
@@ -99,12 +96,9 @@ StatusCode Xgemm<T>::DoGemm(const Layout layout,
   //    matrix A cannot be less than K when rotated, or less than M when not-rotated
   //    matrix B cannot be less than N when rotated, or less than K when not-rotated
   //    matrix C cannot be less than N when rotated, or less than M when not-rotated
-  auto status = TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixB(b_one, b_two, b_buffer, b_offset, b_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixC(c_one, c_two, c_buffer, c_offset, c_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
+  TestMatrixB(b_one, b_two, b_buffer, b_offset, b_ld);
+  TestMatrixC(c_one, c_two, c_buffer, c_offset, c_ld);
 
   // Selects which version of GEMM to run
   const auto do_gemm_direct = (m * n * k < db_["XGEMM_MIN_INDIRECT_SIZE"]);
@@ -131,7 +125,7 @@ StatusCode Xgemm<T>::DoGemm(const Layout layout,
 // requirements, but several pre and post-processing kernels take care of those. However, the
 // overhead of these extra kernels might not be ideal for certain devices/arguments.
 template <typename T>
-StatusCode Xgemm<T>::GemmIndirect(const size_t m, const size_t n, const size_t k,
+void Xgemm<T>::GemmIndirect(const size_t m, const size_t n, const size_t k,
                                   const T alpha,
                                   const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
                                   const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
@@ -142,8 +136,6 @@ StatusCode Xgemm<T>::GemmIndirect(const size_t m, const size_t n, const size_t k
                                   const size_t a_one, const size_t a_two, const bool a_want_rotated,
                                   const size_t b_one, const size_t b_two, const bool b_want_rotated,
                                   const size_t c_one, const size_t c_two, const bool c_want_rotated) {
-  auto status = StatusCode::kSuccess;
-
   // Calculates the ceiled versions of m, n, and k
   const auto m_ceiled = Ceil(m, db_["MWG"]);
   const auto n_ceiled = Ceil(n, db_["NWG"]);
@@ -158,109 +150,95 @@ StatusCode Xgemm<T>::GemmIndirect(const size_t m, const size_t n, const size_t k
   const auto c_one_i = (c_want_rotated) ? n_ceiled : m_ceiled;
   const auto c_two_i = (c_want_rotated) ? m_ceiled : n_ceiled;
 
-  // The padded/transposed input/output matrices: if memory allocation fails, throw an exception
-  try {
-
-    // Loads the program from the database
-    const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-
-    // Determines whether or not temporary matrices are needed
-    auto a_no_temp = a_one == a_one_i && a_two == a_two_i && a_ld == a_one && a_offset == 0 &&
-                     a_do_transpose == false && a_conjugate == false;
-    auto b_no_temp = b_one == b_one_i && b_two == b_two_i && b_ld == b_one && b_offset == 0 &&
-                     b_do_transpose == false && b_conjugate == false;
-    auto c_no_temp = c_one == c_one_i && c_two == c_two_i && c_ld == c_one && c_offset == 0 &&
-                     c_do_transpose == false;
-
-    // Creates the temporary matrices
-    const auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, a_one_i*a_two_i);
-    const auto b_temp = (b_no_temp) ? b_buffer : Buffer<T>(context_, b_one_i*b_two_i);
-    const auto c_temp = (c_no_temp) ? c_buffer : Buffer<T>(context_, c_one_i*c_two_i);
-
-    // Events of all kernels (including pre/post processing kernels)
-    auto eventWaitList = std::vector<Event>();
-    auto emptyEventList = std::vector<Event>();
-
-    // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
-    // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
-    // case nothing has to be done, these kernels can be skipped.
-    if (!a_no_temp) {
-      auto eventProcessA = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
-                                      a_one, a_two, a_ld, a_offset, a_buffer,
-                                      a_one_i, a_two_i, a_one_i, 0, a_temp,
-                                      ConstantOne<T>(), program,
-                                      true, a_do_transpose, a_conjugate);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessA);
-    }
-
-    // As above, but now for matrix B
-    if (!b_no_temp) {
-      auto eventProcessB = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
-                                      b_one, b_two, b_ld, b_offset, b_buffer,
-                                      b_one_i, b_two_i, b_one_i, 0, b_temp,
-                                      ConstantOne<T>(), program,
-                                      true, b_do_transpose, b_conjugate);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessB);
-    }
-
-    // As above, but now for matrix C. This is only necessary if C is used both as input and output.
-    if (!c_no_temp && beta != static_cast<T>(0)) {
-      auto eventProcessC = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
-                                      c_one, c_two, c_ld, c_offset, c_buffer,
-                                      c_one_i, c_two_i, c_one_i, 0, c_temp,
-                                      ConstantOne<T>(), program,
-                                      true, c_do_transpose, false);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessC);
-    }
-
-    // Retrieves the Xgemm kernel from the compiled binary
-    try {
-      auto kernel = Kernel(program, "Xgemm");
-
-      // Sets the kernel arguments
-      kernel.SetArgument(0, static_cast<int>(m_ceiled));
-      kernel.SetArgument(1, static_cast<int>(n_ceiled));
-      kernel.SetArgument(2, static_cast<int>(k_ceiled));
-      kernel.SetArgument(3, GetRealArg(alpha));
-      kernel.SetArgument(4, GetRealArg(beta));
-      kernel.SetArgument(5, a_temp());
-      kernel.SetArgument(6, b_temp());
-      kernel.SetArgument(7, c_temp());
-
-      // Computes the global and local thread sizes
-      const auto global = std::vector<size_t>{
-        (c_one_i * db_["MDIMC"]) / db_["MWG"],
-        (c_two_i * db_["NDIMC"]) / db_["NWG"]
-      };
-      const auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
-
-      // Launches the kernel
-      auto eventKernel = Event();
-      auto eventPointer = (!c_no_temp) ? eventKernel.pointer() : event_;
-      status = RunKernel(kernel, queue_, device_, global, local, eventPointer, eventWaitList);
-      if (ErrorIn(status)) { return status; }
-
-      // Runs the post-processing kernel if needed
-      if (!c_no_temp) {
-        eventWaitList.push_back(eventKernel);
-        status = PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
-                                        c_one_i, c_two_i, c_one_i, 0, c_temp,
-                                        c_one, c_two, c_ld, c_offset, c_buffer,
-                                        ConstantOne<T>(), program,
-                                        false, c_do_transpose, false);
-        if (ErrorIn(status)) { return status; }
-      }
-
-      // Successfully finished the computation
-      return StatusCode::kSuccess;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+  // Loads the program from the database
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+
+  // Determines whether or not temporary matrices are needed
+  auto a_no_temp = a_one == a_one_i && a_two == a_two_i && a_ld == a_one && a_offset == 0 &&
+                   a_do_transpose == false && a_conjugate == false;
+  auto b_no_temp = b_one == b_one_i && b_two == b_two_i && b_ld == b_one && b_offset == 0 &&
+                   b_do_transpose == false && b_conjugate == false;
+  auto c_no_temp = c_one == c_one_i && c_two == c_two_i && c_ld == c_one && c_offset == 0 &&
+                   c_do_transpose == false;
+
+  // Creates the temporary matrices
+  const auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, a_one_i*a_two_i);
+  const auto b_temp = (b_no_temp) ? b_buffer : Buffer<T>(context_, b_one_i*b_two_i);
+  const auto c_temp = (c_no_temp) ? c_buffer : Buffer<T>(context_, c_one_i*c_two_i);
+
+  // Events of all kernels (including pre/post processing kernels)
+  auto eventWaitList = std::vector<Event>();
+  auto emptyEventList = std::vector<Event>();
+
+  // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
+  // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
+  // case nothing has to be done, these kernels can be skipped.
+  if (!a_no_temp) {
+    auto eventProcessA = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
+                           a_one, a_two, a_ld, a_offset, a_buffer,
+                           a_one_i, a_two_i, a_one_i, 0, a_temp,
+                           ConstantOne<T>(), program,
+                           true, a_do_transpose, a_conjugate);
+    eventWaitList.push_back(eventProcessA);
+  }
+
+  // As above, but now for matrix B
+  if (!b_no_temp) {
+    auto eventProcessB = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
+                           b_one, b_two, b_ld, b_offset, b_buffer,
+                           b_one_i, b_two_i, b_one_i, 0, b_temp,
+                           ConstantOne<T>(), program,
+                           true, b_do_transpose, b_conjugate);
+    eventWaitList.push_back(eventProcessB);
+  }
+
+  // As above, but now for matrix C. This is only necessary if C is used both as input and output.
+  if (!c_no_temp && beta != static_cast<T>(0)) {
+    auto eventProcessC = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
+                           c_one, c_two, c_ld, c_offset, c_buffer,
+                           c_one_i, c_two_i, c_one_i, 0, c_temp,
+                           ConstantOne<T>(), program,
+                           true, c_do_transpose, false);
+    eventWaitList.push_back(eventProcessC);
+  }
+
+  // Retrieves the Xgemm kernel from the compiled binary
+  auto kernel = Kernel(program, "Xgemm");
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(m_ceiled));
+  kernel.SetArgument(1, static_cast<int>(n_ceiled));
+  kernel.SetArgument(2, static_cast<int>(k_ceiled));
+  kernel.SetArgument(3, GetRealArg(alpha));
+  kernel.SetArgument(4, GetRealArg(beta));
+  kernel.SetArgument(5, a_temp());
+  kernel.SetArgument(6, b_temp());
+  kernel.SetArgument(7, c_temp());
+
+  // Computes the global and local thread sizes
+  const auto global = std::vector<size_t>{
+    (c_one_i * db_["MDIMC"]) / db_["MWG"],
+    (c_two_i * db_["NDIMC"]) / db_["NWG"]
+  };
+  const auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
+
+  // Launches the kernel
+  auto eventKernel = Event();
+  auto eventPointer = (!c_no_temp) ? eventKernel.pointer() : event_;
+  RunKernel(kernel, queue_, device_, global, local, eventPointer, eventWaitList);
+
+  // Runs the post-processing kernel if needed
+  if (!c_no_temp) {
+    eventWaitList.push_back(eventKernel);
+    PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
+                           c_one_i, c_two_i, c_one_i, 0, c_temp,
+                           c_one, c_two, c_ld, c_offset, c_buffer,
+                           ConstantOne<T>(), program,
+                           false, c_do_transpose, false);
+  }
 }
 
 
@@ -268,7 +246,7 @@ StatusCode Xgemm<T>::GemmIndirect(const size_t m, const size_t n, const size_t k
 
 // The direct version of GEMM, requiring just one kernel, no pre or post-processing kernels.
 template <typename T>
-StatusCode Xgemm<T>::GemmDirect(const size_t m, const size_t n, const size_t k,
+void Xgemm<T>::GemmDirect(const size_t m, const size_t n, const size_t k,
                                 const T alpha,
                                 const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
                                 const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
@@ -281,46 +259,40 @@ StatusCode Xgemm<T>::GemmDirect(const size_t m, const size_t n, const size_t k,
   const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
 
   // Retrieves the proper XgemmDirect kernel from the compiled binary
-  try {
-    const auto name = (a_do_transpose) ? (b_do_transpose ? "XgemmDirectTT" : "XgemmDirectTN") :
-                                         (b_do_transpose ? "XgemmDirectNT" : "XgemmDirectNN");
-    auto kernel = Kernel(program, name);
-
-    // Sets the kernel arguments
-    kernel.SetArgument(0, static_cast<int>(m));
-    kernel.SetArgument(1, static_cast<int>(n));
-    kernel.SetArgument(2, static_cast<int>(k));
-    kernel.SetArgument(3, GetRealArg(alpha));
-    kernel.SetArgument(4, GetRealArg(beta));
-    kernel.SetArgument(5, a_buffer());
-    kernel.SetArgument(6, static_cast<int>(a_offset));
-    kernel.SetArgument(7, static_cast<int>(a_ld));
-    kernel.SetArgument(8, b_buffer());
-    kernel.SetArgument(9, static_cast<int>(b_offset));
-    kernel.SetArgument(10, static_cast<int>(b_ld));
-    kernel.SetArgument(11, c_buffer());
-    kernel.SetArgument(12, static_cast<int>(c_offset));
-    kernel.SetArgument(13, static_cast<int>(c_ld));
-    kernel.SetArgument(14, static_cast<int>(c_do_transpose));
-    kernel.SetArgument(15, static_cast<int>(a_conjugate));
-    kernel.SetArgument(16, static_cast<int>(b_conjugate));
-
-    // Computes the global and local thread sizes
-    const auto m_ceiled = Ceil(m, db_["WGD"]);
-    const auto n_ceiled = Ceil(n, db_["WGD"]);
-    const auto global = std::vector<size_t>{
-      (m_ceiled * db_["MDIMCD"]) / db_["WGD"],
-      (n_ceiled * db_["NDIMCD"]) / db_["WGD"]
-    };
-    const auto local = std::vector<size_t>{db_["MDIMCD"], db_["NDIMCD"]};
-
-    // Launches the kernel
-    auto status = RunKernel(kernel, queue_, device_, global, local, event_);
-    if (ErrorIn(status)) { return status; }
-
-    // Successfully finished the computation
-    return StatusCode::kSuccess;
-  } catch (...) { return StatusCode::kInvalidKernel; }
+  const auto name = (a_do_transpose) ? (b_do_transpose ? "XgemmDirectTT" : "XgemmDirectTN") :
+                                       (b_do_transpose ? "XgemmDirectNT" : "XgemmDirectNN");
+  auto kernel = Kernel(program, name);
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(m));
+  kernel.SetArgument(1, static_cast<int>(n));
+  kernel.SetArgument(2, static_cast<int>(k));
+  kernel.SetArgument(3, GetRealArg(alpha));
+  kernel.SetArgument(4, GetRealArg(beta));
+  kernel.SetArgument(5, a_buffer());
+  kernel.SetArgument(6, static_cast<int>(a_offset));
+  kernel.SetArgument(7, static_cast<int>(a_ld));
+  kernel.SetArgument(8, b_buffer());
+  kernel.SetArgument(9, static_cast<int>(b_offset));
+  kernel.SetArgument(10, static_cast<int>(b_ld));
+  kernel.SetArgument(11, c_buffer());
+  kernel.SetArgument(12, static_cast<int>(c_offset));
+  kernel.SetArgument(13, static_cast<int>(c_ld));
+  kernel.SetArgument(14, static_cast<int>(c_do_transpose));
+  kernel.SetArgument(15, static_cast<int>(a_conjugate));
+  kernel.SetArgument(16, static_cast<int>(b_conjugate));
+
+  // Computes the global and local thread sizes
+  const auto m_ceiled = Ceil(m, db_["WGD"]);
+  const auto n_ceiled = Ceil(n, db_["WGD"]);
+  const auto global = std::vector<size_t>{
+    (m_ceiled * db_["MDIMCD"]) / db_["WGD"],
+    (n_ceiled * db_["NDIMCD"]) / db_["WGD"]
+  };
+  const auto local = std::vector<size_t>{db_["MDIMCD"], db_["NDIMCD"]};
+
+  // Launches the kernel
+  RunKernel(kernel, queue_, device_, global, local, event_);
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xgemm.hpp b/src/routines/level3/xgemm.hpp
index 46e12453..c61611b6 100644
--- a/src/routines/level3/xgemm.hpp
+++ b/src/routines/level3/xgemm.hpp
@@ -28,36 +28,36 @@ class Xgemm: public Routine {
   Xgemm(Queue &queue, EventPointer event, const std::string &name = "GEMM");
 
   // Templated-precision implementation of the routine
-  StatusCode DoGemm(const Layout layout, const Transpose a_transpose, const Transpose b_transpose,
-                    const size_t m, const size_t n, const size_t k,
+  void DoGemm(const Layout layout, const Transpose a_transpose, const Transpose b_transpose,
+              const size_t m, const size_t n, const size_t k,
+              const T alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+              const T beta,
+              const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+
+  // Indirect version of GEMM (with pre and post-processing kernels)
+  void GemmIndirect(const size_t m, const size_t n, const size_t k,
                     const T alpha,
                     const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
                     const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
                     const T beta,
-                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
-
-  // Indirect version of GEMM (with pre and post-processing kernels)
-  StatusCode GemmIndirect(const size_t m, const size_t n, const size_t k,
-                          const T alpha,
-                          const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                          const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                          const T beta,
-                          const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld,
-                          const bool a_do_transpose, const bool b_do_transpose, const bool c_do_transpose,
-                          const bool a_conjugate, const bool b_conjugate,
-                          const size_t a_one, const size_t a_two, const bool a_want_rotated,
-                          const size_t b_one, const size_t b_two, const bool b_want_rotated,
-                          const size_t c_one, const size_t c_two, const bool c_want_rotated);
+                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld,
+                    const bool a_do_transpose, const bool b_do_transpose, const bool c_do_transpose,
+                    const bool a_conjugate, const bool b_conjugate,
+                    const size_t a_one, const size_t a_two, const bool a_want_rotated,
+                    const size_t b_one, const size_t b_two, const bool b_want_rotated,
+                    const size_t c_one, const size_t c_two, const bool c_want_rotated);
 
   // Direct version of GEMM (no pre and post-processing kernels)
-  StatusCode GemmDirect(const size_t m, const size_t n, const size_t k,
-                        const T alpha,
-                        const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                        const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                        const T beta,
-                        const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld,
-                        const bool a_do_transpose, const bool b_do_transpose, const bool c_do_transpose,
-                        const bool a_conjugate, const bool b_conjugate);
+  void GemmDirect(const size_t m, const size_t n, const size_t k,
+                  const T alpha,
+                  const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+                  const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+                  const T beta,
+                  const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld,
+                  const bool a_do_transpose, const bool b_do_transpose, const bool c_do_transpose,
+                  const bool a_conjugate, const bool b_conjugate);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xhemm.cpp b/src/routines/level3/xhemm.cpp
index 9813503e..e5b1502a 100644
--- a/src/routines/level3/xhemm.cpp
+++ b/src/routines/level3/xhemm.cpp
@@ -29,7 +29,7 @@ Xhemm<T>::Xhemm(Queue &queue, EventPointer event, const std::string &name):
 
 // The main routine
 template <typename T>
-StatusCode Xhemm<T>::DoHemm(const Layout layout, const Side side, const Triangle triangle,
+void Xhemm<T>::DoHemm(const Layout layout, const Side side, const Triangle triangle,
                             const size_t m, const size_t n,
                             const T alpha,
                             const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
@@ -38,15 +38,14 @@ StatusCode Xhemm<T>::DoHemm(const Layout layout, const Side side, const Triangle
                             const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((m == 0) || (n == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((m == 0) || (n == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes the k dimension. This is based on whether or not the hermitian matrix is A (on the
   // left) or B (on the right) in the Xgemm routine.
   auto k = (side == Side::kLeft) ? m : n;
 
   // Checks for validity of the squared A matrix
-  auto status = TestMatrixA(k, k, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(k, k, a_buffer, a_offset, a_ld);
 
   // Determines which kernel to run based on the layout (the Xgemm kernel assumes column-major as
   // default) and on whether we are dealing with an upper or lower triangle of the hermitian matrix
@@ -55,73 +54,68 @@ StatusCode Xhemm<T>::DoHemm(const Layout layout, const Side side, const Triangle
   auto kernel_name = (is_upper) ? "HermUpperToSquared" : "HermLowerToSquared";
 
   // Temporary buffer for a copy of the hermitian matrix
-  try {
-    auto temp_herm = Buffer<T>(context_, k*k);
-
-    // Creates a general matrix from the hermitian matrix to be able to run the regular Xgemm
-    // routine afterwards
+  auto temp_herm = Buffer<T>(context_, k*k);
+
+  // Creates a general matrix from the hermitian matrix to be able to run the regular Xgemm
+  // routine afterwards
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the arguments for the hermitian-to-squared kernel
+  kernel.SetArgument(0, static_cast<int>(k));
+  kernel.SetArgument(1, static_cast<int>(a_ld));
+  kernel.SetArgument(2, static_cast<int>(a_offset));
+  kernel.SetArgument(3, a_buffer());
+  kernel.SetArgument(4, static_cast<int>(k));
+  kernel.SetArgument(5, static_cast<int>(k));
+  kernel.SetArgument(6, static_cast<int>(0));
+  kernel.SetArgument(7, temp_herm());
+
+  // Uses the common padding kernel's thread configuration. This is allowed, since the
+  // hermitian-to-squared kernel uses the same parameters.
+  auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
+                                    Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
+  auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
+  auto kernelEvent = Event();
+  RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
+
+  // Synchronize now: 'DoGemm' does not accept a list of events to wait for
+  kernelEvent.WaitForCompletion();
+
+  // Runs the regular Xgemm code with either "C := AB+C" or ...
+  if (side == Side::kLeft) {
+    DoGemm(layout, Transpose::kNo, Transpose::kNo,
+           m, n, k,
+           alpha,
+           temp_herm, 0, k,
+           b_buffer, b_offset, b_ld,
+           beta,
+           c_buffer, c_offset, c_ld);
+  }
+
+  // ... with "C := BA+C". Note that A and B are now reversed.
+  else {
     try {
-      const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the arguments for the hermitian-to-squared kernel
-      kernel.SetArgument(0, static_cast<int>(k));
-      kernel.SetArgument(1, static_cast<int>(a_ld));
-      kernel.SetArgument(2, static_cast<int>(a_offset));
-      kernel.SetArgument(3, a_buffer());
-      kernel.SetArgument(4, static_cast<int>(k));
-      kernel.SetArgument(5, static_cast<int>(k));
-      kernel.SetArgument(6, static_cast<int>(0));
-      kernel.SetArgument(7, temp_herm());
-
-      // Uses the common padding kernel's thread configuration. This is allowed, since the
-      // hermitian-to-squared kernel uses the same parameters.
-      auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
-                                        Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
-      auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
-      auto kernelEvent = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
-      if (ErrorIn(status)) { return status; }
-
-      // Synchronize now: 'DoGemm' does not accept a list of events to wait for
-      kernelEvent.WaitForCompletion();
-
-      // Runs the regular Xgemm code with either "C := AB+C" or ...
-      if (side == Side::kLeft) {
-        status = DoGemm(layout, Transpose::kNo, Transpose::kNo,
-                        m, n, k,
-                        alpha,
-                        temp_herm, 0, k,
-                        b_buffer, b_offset, b_ld,
-                        beta,
-                        c_buffer, c_offset, c_ld);
-      }
-
-      // ... with "C := BA+C". Note that A and B are now reversed.
-      else {
-        status = DoGemm(layout, Transpose::kNo, Transpose::kNo,
-                        m, n, k,
-                        alpha,
-                        b_buffer, b_offset, b_ld,
-                        temp_herm, 0, k,
-                        beta,
-                        c_buffer, c_offset, c_ld);
-
-        // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
-        switch(status) {
-          case StatusCode::kInvalidMatrixA:      status = StatusCode::kInvalidMatrixB; break;
-          case StatusCode::kInvalidMatrixB:      status = StatusCode::kInvalidMatrixA; break;
-          case StatusCode::kInvalidLeadDimA:     status = StatusCode::kInvalidLeadDimB; break;
-          case StatusCode::kInvalidLeadDimB:     status = StatusCode::kInvalidLeadDimA; break;
-          case StatusCode::kInsufficientMemoryA: status = StatusCode::kInsufficientMemoryB; break;
-          case StatusCode::kInsufficientMemoryB: status = StatusCode::kInsufficientMemoryA; break;
-        }
+      DoGemm(layout, Transpose::kNo, Transpose::kNo,
+             m, n, k,
+             alpha,
+             b_buffer, b_offset, b_ld,
+             temp_herm, 0, k,
+             beta,
+             c_buffer, c_offset, c_ld);
+    } catch (BLASError &e) {
+      // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
+      switch(e.status()) {
+        case StatusCode::kInvalidMatrixA:      throw BLASError(StatusCode::kInvalidMatrixB, e.details());
+        case StatusCode::kInvalidMatrixB:      throw BLASError(StatusCode::kInvalidMatrixA, e.details());
+        case StatusCode::kInvalidLeadDimA:     throw BLASError(StatusCode::kInvalidLeadDimB, e.details());
+        case StatusCode::kInvalidLeadDimB:     throw BLASError(StatusCode::kInvalidLeadDimA, e.details());
+        case StatusCode::kInsufficientMemoryA: throw BLASError(StatusCode::kInsufficientMemoryB, e.details());
+        case StatusCode::kInsufficientMemoryB: throw BLASError(StatusCode::kInsufficientMemoryA, e.details());
+        default:                               throw;
       }
-
-      // Return the status of the Xgemm routine
-      return status;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+    }
+  }
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xhemm.hpp b/src/routines/level3/xhemm.hpp
index 272bd2ec..2385706e 100644
--- a/src/routines/level3/xhemm.hpp
+++ b/src/routines/level3/xhemm.hpp
@@ -37,13 +37,13 @@ class Xhemm: public Xgemm<T> {
   Xhemm(Queue &queue, EventPointer event, const std::string &name = "HEMM");
 
   // Templated-precision implementation of the routine
-  StatusCode DoHemm(const Layout layout, const Side side, const Triangle triangle,
-                    const size_t m, const size_t n,
-                    const T alpha,
-                    const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                    const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                    const T beta,
-                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoHemm(const Layout layout, const Side side, const Triangle triangle,
+              const size_t m, const size_t n,
+              const T alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+              const T beta,
+              const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xher2k.cpp b/src/routines/level3/xher2k.cpp
index bf328729..ee3bb8b8 100644
--- a/src/routines/level3/xher2k.cpp
+++ b/src/routines/level3/xher2k.cpp
@@ -22,8 +22,7 @@ namespace clblast {
 // Constructor: forwards to base class constructor
 template <typename T, typename U>
 Xher2k<T,U>::Xher2k(Queue &queue, EventPointer event, const std::string &name):
-    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>()) {
-  source_string_ =
+    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>(), {}, {
     #include "../../kernels/level3/level3.opencl"
     #include "../../kernels/level3/copy_fast.opencl"
     #include "../../kernels/level3/copy_pad.opencl"
@@ -32,23 +31,23 @@ Xher2k<T,U>::Xher2k(Queue &queue, EventPointer event, const std::string &name):
     #include "../../kernels/level3/xgemm_part1.opencl"
     #include "../../kernels/level3/xgemm_part2.opencl"
     #include "../../kernels/level3/xgemm_part3.opencl"
-  ;
+    }) {
 }
 
 // =================================================================================================
 
 // The main routine
 template <typename T, typename U>
-StatusCode Xher2k<T,U>::DoHer2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
-                                const size_t n, const size_t k,
-                                const T alpha,
-                                const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                                const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                                const U beta,
-                                const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
+void Xher2k<T,U>::DoHer2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
+                          const size_t n, const size_t k,
+                          const T alpha,
+                          const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+                          const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+                          const U beta,
+                          const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((n == 0) || (k == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((n == 0) || (k == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Determines whether to apply the conjugate transpose to matrix B (argument: no transpose) or
   // to matrix A (argument: conjugate transpose)
@@ -71,12 +70,9 @@ StatusCode Xher2k<T,U>::DoHer2k(const Layout layout, const Triangle triangle, co
   //    matrix A cannot be less than N when rotated, or less than K when not-rotated
   //    matrix B cannot be less than N when rotated, or less than K when not-rotated
   //    matrix C cannot be less than N
-  auto status = TestMatrixA(ab_one, ab_two, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixB(ab_one, ab_two, b_buffer, b_offset, b_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixC(n, n, c_buffer, c_offset, c_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(ab_one, ab_two, a_buffer, a_offset, a_ld);
+  TestMatrixB(ab_one, ab_two, b_buffer, b_offset, b_ld);
+  TestMatrixC(n, n, c_buffer, c_offset, c_ld);
 
   // Calculates the ceiled versions of n and k
   auto n_ceiled = Ceil(Ceil(n, db_["MWG"]), db_["NWG"]);
@@ -85,145 +81,128 @@ StatusCode Xher2k<T,U>::DoHer2k(const Layout layout, const Triangle triangle, co
   // Decides which kernel to run: the upper-triangular or lower-triangular version
   auto kernel_name = (triangle == Triangle::kUpper) ? "XgemmUpper" : "XgemmLower";
 
-  // The padded/transposed input/output matrices: if memory allocation fails, throw an exception
-  try {
-
-    // Loads the program from the database
-    const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-
-    // Determines whether or not temporary matrices are needed
-    auto a1_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                      ab_rotated == false && ab_conjugate == false;
-    auto a2_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                      ab_rotated == false && ab_conjugate == true;
-    auto b1_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
-                      ab_rotated == false && ab_conjugate == false;
-    auto b2_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
-                      ab_rotated == false && ab_conjugate == true;
-
-    // Creates the temporary matrices
-    auto a1_temp = (a1_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto a2_temp = (a2_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto b1_temp = (b1_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto b2_temp = (b2_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
-
-    // Convert the arguments to complex versions
-    auto complex_beta = T{beta, static_cast<U>(0.0)};
-
-    // Events of all kernels (including pre/post processing kernels)
-    auto eventWaitList = std::vector<Event>();
-    auto emptyEventList = std::vector<Event>();
-
-    // Runs the pre-processing kernels. This transposes the matrices A and B, but also pads zeros to
-    // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
-    // case nothing has to be done, these kernels can be skipped.
-    if (!a1_no_temp) {
-      auto eventProcessA1 = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA1.pointer(), emptyEventList,
-                                      ab_one, ab_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, a1_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, ab_conjugate);
-      eventWaitList.push_back(eventProcessA1);
-      if (ErrorIn(status)) { return status; }
-    }
-    if (!a2_no_temp) {
-      auto eventProcessA2 = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA2.pointer(), emptyEventList,
-                                      ab_one, ab_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, a2_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, !ab_conjugate);
-      eventWaitList.push_back(eventProcessA2);
-      if (ErrorIn(status)) { return status; }
-    }
-    if (!b1_no_temp) {
-      auto eventProcessB1 = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB1.pointer(), emptyEventList,
-                                      ab_one, ab_two, b_ld, b_offset, b_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, b1_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, ab_conjugate);
-      eventWaitList.push_back(eventProcessB1);
-      if (ErrorIn(status)) { return status; }
-    }
-    if (!b2_no_temp) {
-      auto eventProcessB2 = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB2.pointer(), emptyEventList,
-                                      ab_one, ab_two, b_ld, b_offset, b_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, b2_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, !ab_conjugate);
-      eventWaitList.push_back(eventProcessB2);
-      if (ErrorIn(status)) { return status; }
-    }
-
-    // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
-    // modify the other triangle.
-    auto eventProcessC = Event();
-    status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
-                                    n, n, c_ld, c_offset, c_buffer,
-                                    n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                    ConstantOne<T>(), program,
-                                    true, c_rotated, false);
-    eventWaitList.push_back(eventProcessC);
-    if (ErrorIn(status)) { return status; }
-
-    // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
-    try {
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the kernel arguments
-      kernel.SetArgument(0, static_cast<int>(n_ceiled));
-      kernel.SetArgument(1, static_cast<int>(k_ceiled));
-      kernel.SetArgument(2, GetRealArg(alpha));
-      kernel.SetArgument(3, GetRealArg(complex_beta));
-      kernel.SetArgument(4, a1_temp());
-      kernel.SetArgument(5, b2_temp());
-      kernel.SetArgument(6, c_temp());
-
-      // Computes the global and local thread sizes
-      auto global = std::vector<size_t>{
-        (n_ceiled * db_["MDIMC"]) / db_["MWG"],
-        (n_ceiled * db_["NDIMC"]) / db_["NWG"]
-      };
-      auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
-
-      // Launches the kernel
-      auto eventKernel1 = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel1.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel1);
-
-      // Swaps the arguments for matrices A and B, sets 'beta' to 1, and conjugate alpha
-      auto conjugate_alpha = T{alpha.real(), -alpha.imag()};
-      auto complex_one = T{static_cast<U>(1.0), static_cast<U>(0.0)};
-      kernel.SetArgument(2, GetRealArg(conjugate_alpha));
-      kernel.SetArgument(3, GetRealArg(complex_one));
-      kernel.SetArgument(4, b1_temp());
-      kernel.SetArgument(5, a2_temp());
-
-      // Runs the kernel again
-      auto eventKernel2 = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel2.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel2);
-
-      // Runs the post-processing kernel
-      auto upper = (triangle == Triangle::kUpper);
-      auto lower = (triangle == Triangle::kLower);
-      status = PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
-                                      n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                      n, n, c_ld, c_offset, c_buffer,
-                                      ConstantOne<T>(), program,
-                                      false, c_rotated, false, upper, lower, true);
-      if (ErrorIn(status)) { return status; }
-
-      // Successfully finished the computation
-      return StatusCode::kSuccess;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+  // Loads the program from the database
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+
+  // Determines whether or not temporary matrices are needed
+  auto a1_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                    ab_rotated == false && ab_conjugate == false;
+  auto a2_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                    ab_rotated == false && ab_conjugate == true;
+  auto b1_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
+                    ab_rotated == false && ab_conjugate == false;
+  auto b2_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
+                    ab_rotated == false && ab_conjugate == true;
+
+  // Creates the temporary matrices
+  auto a1_temp = (a1_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto a2_temp = (a2_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto b1_temp = (b1_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto b2_temp = (b2_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
+
+  // Convert the arguments to complex versions
+  auto complex_beta = T{beta, static_cast<U>(0.0)};
+
+  // Events of all kernels (including pre/post processing kernels)
+  auto eventWaitList = std::vector<Event>();
+  auto emptyEventList = std::vector<Event>();
+
+  // Runs the pre-processing kernels. This transposes the matrices A and B, but also pads zeros to
+  // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
+  // case nothing has to be done, these kernels can be skipped.
+  if (!a1_no_temp) {
+    auto eventProcessA1 = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA1.pointer(), emptyEventList,
+                           ab_one, ab_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, a1_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, ab_conjugate);
+    eventWaitList.push_back(eventProcessA1);
+  }
+  if (!a2_no_temp) {
+    auto eventProcessA2 = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA2.pointer(), emptyEventList,
+                           ab_one, ab_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, a2_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, !ab_conjugate);
+    eventWaitList.push_back(eventProcessA2);
+  }
+  if (!b1_no_temp) {
+    auto eventProcessB1 = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB1.pointer(), emptyEventList,
+                           ab_one, ab_two, b_ld, b_offset, b_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, b1_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, ab_conjugate);
+    eventWaitList.push_back(eventProcessB1);
+  }
+  if (!b2_no_temp) {
+    auto eventProcessB2 = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB2.pointer(), emptyEventList,
+                           ab_one, ab_two, b_ld, b_offset, b_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, b2_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, !ab_conjugate);
+    eventWaitList.push_back(eventProcessB2);
+  }
+
+  // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
+  // modify the other triangle.
+  auto eventProcessC = Event();
+  PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
+                         n, n, c_ld, c_offset, c_buffer,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         ConstantOne<T>(), program,
+                         true, c_rotated, false);
+  eventWaitList.push_back(eventProcessC);
+
+  // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(n_ceiled));
+  kernel.SetArgument(1, static_cast<int>(k_ceiled));
+  kernel.SetArgument(2, GetRealArg(alpha));
+  kernel.SetArgument(3, GetRealArg(complex_beta));
+  kernel.SetArgument(4, a1_temp());
+  kernel.SetArgument(5, b2_temp());
+  kernel.SetArgument(6, c_temp());
+
+  // Computes the global and local thread sizes
+  auto global = std::vector<size_t>{
+    (n_ceiled * db_["MDIMC"]) / db_["MWG"],
+    (n_ceiled * db_["NDIMC"]) / db_["NWG"]
+  };
+  auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
+
+  // Launches the kernel
+  auto eventKernel1 = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel1.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel1);
+
+  // Swaps the arguments for matrices A and B, sets 'beta' to 1, and conjugate alpha
+  auto conjugate_alpha = T{alpha.real(), -alpha.imag()};
+  auto complex_one = T{static_cast<U>(1.0), static_cast<U>(0.0)};
+  kernel.SetArgument(2, GetRealArg(conjugate_alpha));
+  kernel.SetArgument(3, GetRealArg(complex_one));
+  kernel.SetArgument(4, b1_temp());
+  kernel.SetArgument(5, a2_temp());
+
+  // Runs the kernel again
+  auto eventKernel2 = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel2.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel2);
+
+  // Runs the post-processing kernel
+  auto upper = (triangle == Triangle::kUpper);
+  auto lower = (triangle == Triangle::kLower);
+  PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         n, n, c_ld, c_offset, c_buffer,
+                         ConstantOne<T>(), program,
+                         false, c_rotated, false, upper, lower, true);
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xher2k.hpp b/src/routines/level3/xher2k.hpp
index 23996219..acc346e4 100644
--- a/src/routines/level3/xher2k.hpp
+++ b/src/routines/level3/xher2k.hpp
@@ -30,13 +30,13 @@ class Xher2k: public Routine {
   Xher2k(Queue &queue, EventPointer event, const std::string &name = "HER2K");
 
   // Templated-precision implementation of the routine
-  StatusCode DoHer2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
-                     const size_t n, const size_t k,
-                     const T alpha,
-                     const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                     const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                     const U beta,
-                     const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoHer2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
+               const size_t n, const size_t k,
+               const T alpha,
+               const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+               const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+               const U beta,
+               const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xherk.cpp b/src/routines/level3/xherk.cpp
index 77422526..ae8e9324 100644
--- a/src/routines/level3/xherk.cpp
+++ b/src/routines/level3/xherk.cpp
@@ -22,8 +22,7 @@ namespace clblast {
 // Constructor: forwards to base class constructor
 template <typename T, typename U>
 Xherk<T,U>::Xherk(Queue &queue, EventPointer event, const std::string &name):
-    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>()) {
-  source_string_ =
+    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>(), {}, {
     #include "../../kernels/level3/level3.opencl"
     #include "../../kernels/level3/copy_fast.opencl"
     #include "../../kernels/level3/copy_pad.opencl"
@@ -32,14 +31,14 @@ Xherk<T,U>::Xherk(Queue &queue, EventPointer event, const std::string &name):
     #include "../../kernels/level3/xgemm_part1.opencl"
     #include "../../kernels/level3/xgemm_part2.opencl"
     #include "../../kernels/level3/xgemm_part3.opencl"
-  ;
+    }) {
 }
 
 // =================================================================================================
 
 // The main routine
 template <typename T, typename U>
-StatusCode Xherk<T,U>::DoHerk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
+void Xherk<T,U>::DoHerk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
                               const size_t n, const size_t k,
                               const U alpha,
                               const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
@@ -47,7 +46,7 @@ StatusCode Xherk<T,U>::DoHerk(const Layout layout, const Triangle triangle, cons
                               const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((n == 0) || (k == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((n == 0) || (k == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Determines whether to apply the conjugate transpose to matrix B (argument: no transpose) or
   // to matrix A (argument: conjugate transpose)
@@ -70,10 +69,8 @@ StatusCode Xherk<T,U>::DoHerk(const Layout layout, const Triangle triangle, cons
   // space. Also tests that the leading dimensions of:
   //    matrix A cannot be less than N when rotated, or less than K when not-rotated
   //    matrix C cannot be less than N
-  auto status = TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixC(n, n, c_buffer, c_offset, c_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
+  TestMatrixC(n, n, c_buffer, c_offset, c_ld);
 
   // Calculates the ceiled versions of n and k
   auto n_ceiled = Ceil(Ceil(n, db_["MWG"]), db_["NWG"]);
@@ -82,106 +79,92 @@ StatusCode Xherk<T,U>::DoHerk(const Layout layout, const Triangle triangle, cons
   // Decides which kernel to run: the upper-triangular or lower-triangular version
   auto kernel_name = (triangle == Triangle::kUpper) ? "XgemmUpper" : "XgemmLower";
 
-  // The padded/transposed input/output matrices: if memory allocation fails, throw an exception
-  try {
-
-    // Loads the program from the database
-    const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-
-    // Determines whether or not temporary matrices are needed
-    auto a_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                     a_rotated == false && a_conjugate == false;
-    auto b_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                     a_rotated == false && b_conjugate == false;
-
-    // Creates the temporary matrices
-    auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto b_temp = (b_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
-
-    // Convert the arguments to complex versions
-    auto complex_alpha = T{alpha, static_cast<U>(0.0)};
-    auto complex_beta = T{beta, static_cast<U>(0.0)};
-
-    // Events of all kernels (including pre/post processing kernels)
-    auto eventWaitList = std::vector<Event>();
-    auto emptyEventList = std::vector<Event>();
-
-    // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
-    // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
-    // case nothing has to be done, these kernels can be skipped. Two copies are created.
-    if (!a_no_temp) {
-      auto eventProcessA = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
-                                      a_one, a_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
-                                      ConstantOne<T>(), program,
-                                      true, a_rotated, a_conjugate);
-      eventWaitList.push_back(eventProcessA);
-      if (ErrorIn(status)) { return status; }
-    }
-    if (!b_no_temp) {
-      auto eventProcessB = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
-                                      a_one, a_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, b_temp,
-                                      ConstantOne<T>(), program,
-                                      true, a_rotated, b_conjugate);
-      eventWaitList.push_back(eventProcessB);
-      if (ErrorIn(status)) { return status; }
-    }
-
-    // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
-    // modify the other triangle.
-    auto eventProcessC = Event();
-    status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
-                                    n, n, c_ld, c_offset, c_buffer,
-                                    n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                    ConstantOne<T>(), program,
-                                    true, c_rotated, false);
-    eventWaitList.push_back(eventProcessC);
-    if (ErrorIn(status)) { return status; }
-
-    // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
-    try {
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the kernel arguments
-      kernel.SetArgument(0, static_cast<int>(n_ceiled));
-      kernel.SetArgument(1, static_cast<int>(k_ceiled));
-      kernel.SetArgument(2, GetRealArg(complex_alpha));
-      kernel.SetArgument(3, GetRealArg(complex_beta));
-      kernel.SetArgument(4, a_temp());
-      kernel.SetArgument(5, b_temp());
-      kernel.SetArgument(6, c_temp());
-
-      // Computes the global and local thread sizes
-      auto global = std::vector<size_t>{
-        (n_ceiled * db_["MDIMC"]) / db_["MWG"],
-        (n_ceiled * db_["NDIMC"]) / db_["NWG"]
-      };
-      auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
-
-      // Launches the kernel
-      auto eventKernel = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel);
-
-      // Runs the post-processing kernel
-      auto upper = (triangle == Triangle::kUpper);
-      auto lower = (triangle == Triangle::kLower);
-      status = PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
-                                      n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                      n, n, c_ld, c_offset, c_buffer,
-                                      ConstantOne<T>(), program,
-                                      false, c_rotated, false, upper, lower, true);
-      if (ErrorIn(status)) { return status; }
-
-      // Successfully finished the computation
-      return StatusCode::kSuccess;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+  // Loads the program from the database
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+
+  // Determines whether or not temporary matrices are needed
+  auto a_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                   a_rotated == false && a_conjugate == false;
+  auto b_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                   a_rotated == false && b_conjugate == false;
+
+  // Creates the temporary matrices
+  auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto b_temp = (b_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
+
+  // Convert the arguments to complex versions
+  auto complex_alpha = T{alpha, static_cast<U>(0.0)};
+  auto complex_beta = T{beta, static_cast<U>(0.0)};
+
+  // Events of all kernels (including pre/post processing kernels)
+  auto eventWaitList = std::vector<Event>();
+  auto emptyEventList = std::vector<Event>();
+
+  // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
+  // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
+  // case nothing has to be done, these kernels can be skipped. Two copies are created.
+  if (!a_no_temp) {
+    auto eventProcessA = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
+                           a_one, a_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
+                           ConstantOne<T>(), program,
+                           true, a_rotated, a_conjugate);
+    eventWaitList.push_back(eventProcessA);
+  }
+  if (!b_no_temp) {
+    auto eventProcessB = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
+                           a_one, a_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, b_temp,
+                           ConstantOne<T>(), program,
+                           true, a_rotated, b_conjugate);
+    eventWaitList.push_back(eventProcessB);
+  }
+
+  // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
+  // modify the other triangle.
+  auto eventProcessC = Event();
+  PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
+                         n, n, c_ld, c_offset, c_buffer,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         ConstantOne<T>(), program,
+                         true, c_rotated, false);
+  eventWaitList.push_back(eventProcessC);
+
+  // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(n_ceiled));
+  kernel.SetArgument(1, static_cast<int>(k_ceiled));
+  kernel.SetArgument(2, GetRealArg(complex_alpha));
+  kernel.SetArgument(3, GetRealArg(complex_beta));
+  kernel.SetArgument(4, a_temp());
+  kernel.SetArgument(5, b_temp());
+  kernel.SetArgument(6, c_temp());
+
+  // Computes the global and local thread sizes
+  auto global = std::vector<size_t>{
+    (n_ceiled * db_["MDIMC"]) / db_["MWG"],
+    (n_ceiled * db_["NDIMC"]) / db_["NWG"]
+  };
+  auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
+
+  // Launches the kernel
+  auto eventKernel = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel);
+
+  // Runs the post-processing kernel
+  auto upper = (triangle == Triangle::kUpper);
+  auto lower = (triangle == Triangle::kLower);
+  PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         n, n, c_ld, c_offset, c_buffer,
+                         ConstantOne<T>(), program,
+                         false, c_rotated, false, upper, lower, true);
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xherk.hpp b/src/routines/level3/xherk.hpp
index 3f156a1b..51f29d7e 100644
--- a/src/routines/level3/xherk.hpp
+++ b/src/routines/level3/xherk.hpp
@@ -30,12 +30,12 @@ class Xherk: public Routine {
   Xherk(Queue &queue, EventPointer event, const std::string &name = "HERK");
 
   // Templated-precision implementation of the routine
-  StatusCode DoHerk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
-                    const size_t n, const size_t k,
-                    const U alpha,
-                    const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                    const U beta,
-                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoHerk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
+              const size_t n, const size_t k,
+              const U alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const U beta,
+              const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xsymm.cpp b/src/routines/level3/xsymm.cpp
index 04e4b718..d7f771d1 100644
--- a/src/routines/level3/xsymm.cpp
+++ b/src/routines/level3/xsymm.cpp
@@ -29,7 +29,7 @@ Xsymm<T>::Xsymm(Queue &queue, EventPointer event, const std::string &name):
 
 // The main routine
 template <typename T>
-StatusCode Xsymm<T>::DoSymm(const Layout layout, const Side side, const Triangle triangle,
+void Xsymm<T>::DoSymm(const Layout layout, const Side side, const Triangle triangle,
                             const size_t m, const size_t n,
                             const T alpha,
                             const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
@@ -38,15 +38,14 @@ StatusCode Xsymm<T>::DoSymm(const Layout layout, const Side side, const Triangle
                             const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((m == 0) || (n == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((m == 0) || (n == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes the k dimension. This is based on whether or not the symmetric matrix is A (on the
   // left) or B (on the right) in the Xgemm routine.
   auto k = (side == Side::kLeft) ? m : n;
 
   // Checks for validity of the squared A matrix
-  auto status = TestMatrixA(k, k, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(k, k, a_buffer, a_offset, a_ld);
 
   // Determines which kernel to run based on the layout (the Xgemm kernel assumes column-major as
   // default) and on whether we are dealing with an upper or lower triangle of the symmetric matrix
@@ -55,73 +54,68 @@ StatusCode Xsymm<T>::DoSymm(const Layout layout, const Side side, const Triangle
   auto kernel_name = (is_upper) ? "SymmUpperToSquared" : "SymmLowerToSquared";
 
   // Temporary buffer for a copy of the symmetric matrix
-  try {
-    auto temp_symm = Buffer<T>(context_, k*k);
-
-    // Creates a general matrix from the symmetric matrix to be able to run the regular Xgemm
-    // routine afterwards
+  auto temp_symm = Buffer<T>(context_, k*k);
+
+  // Creates a general matrix from the symmetric matrix to be able to run the regular Xgemm
+  // routine afterwards
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the arguments for the symmetric-to-squared kernel
+  kernel.SetArgument(0, static_cast<int>(k));
+  kernel.SetArgument(1, static_cast<int>(a_ld));
+  kernel.SetArgument(2, static_cast<int>(a_offset));
+  kernel.SetArgument(3, a_buffer());
+  kernel.SetArgument(4, static_cast<int>(k));
+  kernel.SetArgument(5, static_cast<int>(k));
+  kernel.SetArgument(6, static_cast<int>(0));
+  kernel.SetArgument(7, temp_symm());
+
+  // Uses the common padding kernel's thread configuration. This is allowed, since the
+  // symmetric-to-squared kernel uses the same parameters.
+  auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
+                                    Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
+  auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
+  auto kernelEvent = Event();
+  RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
+
+  // Synchronize now: 'DoGemm' does not accept a list of events to wait for
+  kernelEvent.WaitForCompletion();
+
+  // Runs the regular Xgemm code with either "C := AB+C" or ...
+  if (side == Side::kLeft) {
+    DoGemm(layout, Transpose::kNo, Transpose::kNo,
+           m, n, k,
+           alpha,
+           temp_symm, 0, k,
+           b_buffer, b_offset, b_ld,
+           beta,
+           c_buffer, c_offset, c_ld);
+  }
+
+  // ... with "C := BA+C". Note that A and B are now reversed.
+  else {
     try {
-      const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the arguments for the symmetric-to-squared kernel
-      kernel.SetArgument(0, static_cast<int>(k));
-      kernel.SetArgument(1, static_cast<int>(a_ld));
-      kernel.SetArgument(2, static_cast<int>(a_offset));
-      kernel.SetArgument(3, a_buffer());
-      kernel.SetArgument(4, static_cast<int>(k));
-      kernel.SetArgument(5, static_cast<int>(k));
-      kernel.SetArgument(6, static_cast<int>(0));
-      kernel.SetArgument(7, temp_symm());
-
-      // Uses the common padding kernel's thread configuration. This is allowed, since the
-      // symmetric-to-squared kernel uses the same parameters.
-      auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
-                                        Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
-      auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
-      auto kernelEvent = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
-      if (ErrorIn(status)) { return status; }
-
-      // Synchronize now: 'DoGemm' does not accept a list of events to wait for
-      kernelEvent.WaitForCompletion();
-
-      // Runs the regular Xgemm code with either "C := AB+C" or ...
-      if (side == Side::kLeft) {
-        status = DoGemm(layout, Transpose::kNo, Transpose::kNo,
-                        m, n, k,
-                        alpha,
-                        temp_symm, 0, k,
-                        b_buffer, b_offset, b_ld,
-                        beta,
-                        c_buffer, c_offset, c_ld);
-      }
-
-      // ... with "C := BA+C". Note that A and B are now reversed.
-      else {
-        status = DoGemm(layout, Transpose::kNo, Transpose::kNo,
-                        m, n, k,
-                        alpha,
-                        b_buffer, b_offset, b_ld,
-                        temp_symm, 0, k,
-                        beta,
-                        c_buffer, c_offset, c_ld);
-
-        // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
-        switch(status) {
-          case StatusCode::kInvalidMatrixA:      status = StatusCode::kInvalidMatrixB; break;
-          case StatusCode::kInvalidMatrixB:      status = StatusCode::kInvalidMatrixA; break;
-          case StatusCode::kInvalidLeadDimA:     status = StatusCode::kInvalidLeadDimB; break;
-          case StatusCode::kInvalidLeadDimB:     status = StatusCode::kInvalidLeadDimA; break;
-          case StatusCode::kInsufficientMemoryA: status = StatusCode::kInsufficientMemoryB; break;
-          case StatusCode::kInsufficientMemoryB: status = StatusCode::kInsufficientMemoryA; break;
-        }
+      DoGemm(layout, Transpose::kNo, Transpose::kNo,
+             m, n, k,
+             alpha,
+             b_buffer, b_offset, b_ld,
+             temp_symm, 0, k,
+             beta,
+             c_buffer, c_offset, c_ld);
+    } catch (BLASError &e) {
+      // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
+      switch(e.status()) {
+        case StatusCode::kInvalidMatrixA:      throw BLASError(StatusCode::kInvalidMatrixB, e.details());
+        case StatusCode::kInvalidMatrixB:      throw BLASError(StatusCode::kInvalidMatrixA, e.details());
+        case StatusCode::kInvalidLeadDimA:     throw BLASError(StatusCode::kInvalidLeadDimB, e.details());
+        case StatusCode::kInvalidLeadDimB:     throw BLASError(StatusCode::kInvalidLeadDimA, e.details());
+        case StatusCode::kInsufficientMemoryA: throw BLASError(StatusCode::kInsufficientMemoryB, e.details());
+        case StatusCode::kInsufficientMemoryB: throw BLASError(StatusCode::kInsufficientMemoryA, e.details());
+        default:                               throw;
       }
-
-      // Return the status of the Xgemm routine
-      return status;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+    }
+  }
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xsymm.hpp b/src/routines/level3/xsymm.hpp
index 428f78ef..ee965364 100644
--- a/src/routines/level3/xsymm.hpp
+++ b/src/routines/level3/xsymm.hpp
@@ -39,13 +39,13 @@ class Xsymm: public Xgemm<T> {
   Xsymm(Queue &queue, EventPointer event, const std::string &name = "SYMM");
 
   // Templated-precision implementation of the routine
-  StatusCode DoSymm(const Layout layout, const Side side, const Triangle triangle,
-                    const size_t m, const size_t n,
-                    const T alpha,
-                    const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                    const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                    const T beta,
-                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoSymm(const Layout layout, const Side side, const Triangle triangle,
+              const size_t m, const size_t n,
+              const T alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+              const T beta,
+              const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xsyr2k.cpp b/src/routines/level3/xsyr2k.cpp
index badf3100..cb0e0461 100644
--- a/src/routines/level3/xsyr2k.cpp
+++ b/src/routines/level3/xsyr2k.cpp
@@ -22,8 +22,7 @@ namespace clblast {
 // Constructor: forwards to base class constructor
 template <typename T>
 Xsyr2k<T>::Xsyr2k(Queue &queue, EventPointer event, const std::string &name):
-    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>()) {
-  source_string_ =
+    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>(), {}, {
     #include "../../kernels/level3/level3.opencl"
     #include "../../kernels/level3/copy_fast.opencl"
     #include "../../kernels/level3/copy_pad.opencl"
@@ -32,14 +31,14 @@ Xsyr2k<T>::Xsyr2k(Queue &queue, EventPointer event, const std::string &name):
     #include "../../kernels/level3/xgemm_part1.opencl"
     #include "../../kernels/level3/xgemm_part2.opencl"
     #include "../../kernels/level3/xgemm_part3.opencl"
-  ;
+    }) {
 }
 
 // =================================================================================================
 
 // The main routine
 template <typename T>
-StatusCode Xsyr2k<T>::DoSyr2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
+void Xsyr2k<T>::DoSyr2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
                               const size_t n, const size_t k,
                               const T alpha,
                               const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
@@ -48,7 +47,7 @@ StatusCode Xsyr2k<T>::DoSyr2k(const Layout layout, const Triangle triangle, cons
                               const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((n == 0) || (k == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((n == 0) || (k == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes whether or not the matrices are transposed in memory. This is based on their layout
   // (row or column-major) and whether or not they are requested to be pre-transposed.
@@ -67,12 +66,9 @@ StatusCode Xsyr2k<T>::DoSyr2k(const Layout layout, const Triangle triangle, cons
   //    matrix A cannot be less than N when rotated, or less than K when not-rotated
   //    matrix B cannot be less than N when rotated, or less than K when not-rotated
   //    matrix C cannot be less than N
-  auto status = TestMatrixA(ab_one, ab_two, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixB(ab_one, ab_two, b_buffer, b_offset, b_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixC(n, n, c_buffer, c_offset, c_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(ab_one, ab_two, a_buffer, a_offset, a_ld);
+  TestMatrixB(ab_one, ab_two, b_buffer, b_offset, b_ld);
+  TestMatrixC(n, n, c_buffer, c_offset, c_ld);
 
   // Calculates the ceiled versions of n and k
   auto n_ceiled = Ceil(Ceil(n, db_["MWG"]), db_["NWG"]);
@@ -81,114 +77,99 @@ StatusCode Xsyr2k<T>::DoSyr2k(const Layout layout, const Triangle triangle, cons
   // Decides which kernel to run: the upper-triangular or lower-triangular version
   auto kernel_name = (triangle == Triangle::kUpper) ? "XgemmUpper" : "XgemmLower";
 
-  // The padded/transposed input/output matrices: if memory allocation fails, throw an exception
-  try {
-
-    // Loads the program from the database
-    const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-
-    // Determines whether or not temporary matrices are needed
-    auto a_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                     ab_rotated == false;
-    auto b_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
-                     ab_rotated == false;
-
-    // Creates the temporary matrices
-    auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto b_temp = (b_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
-
-    // Events of all kernels (including pre/post processing kernels)
-    auto eventWaitList = std::vector<Event>();
-    auto emptyEventList = std::vector<Event>();
-
-    // Runs the pre-processing kernels. This transposes the matrices A and B, but also pads zeros to
-    // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
-    // case nothing has to be done, these kernels can be skipped.
-    if (!a_no_temp) {
-      auto eventProcessA = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
-                                      ab_one, ab_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, false);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessA);
-    }
-    if (!b_no_temp) {
-      auto eventProcessB = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
-                                      ab_one, ab_two, b_ld, b_offset, b_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, b_temp,
-                                      ConstantOne<T>(), program,
-                                      true, ab_rotated, false);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessB);
-    }
-
-    // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
-    // modify the other triangle.
-    auto eventProcessC = Event();
-    status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
-                                    n, n, c_ld, c_offset, c_buffer,
-                                    n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                    ConstantOne<T>(), program,
-                                    true, c_rotated, false);
-    if (ErrorIn(status)) { return status; }
-    eventWaitList.push_back(eventProcessC);
-
-    // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
-    try {
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the kernel arguments
-      kernel.SetArgument(0, static_cast<int>(n_ceiled));
-      kernel.SetArgument(1, static_cast<int>(k_ceiled));
-      kernel.SetArgument(2, GetRealArg(alpha));
-      kernel.SetArgument(3, GetRealArg(beta));
-      kernel.SetArgument(4, a_temp());
-      kernel.SetArgument(5, b_temp());
-      kernel.SetArgument(6, c_temp());
-
-      // Computes the global and local thread sizes
-      auto global = std::vector<size_t>{
-        (n_ceiled * db_["MDIMC"]) / db_["MWG"],
-        (n_ceiled * db_["NDIMC"]) / db_["NWG"]
-      };
-      auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
-
-      // Launches the kernel
-      auto eventKernel1 = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel1.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel1);
-
-      // Swaps the arguments for matrices A and B, and sets 'beta' to 1
-      auto one = static_cast<T>(1);
-      kernel.SetArgument(3, GetRealArg(one));
-      kernel.SetArgument(4, b_temp());
-      kernel.SetArgument(5, a_temp());
-
-      // Runs the kernel again
-      auto eventKernel2 = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel2.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel2);
-
-      // Runs the post-processing kernel
-      auto upper = (triangle == Triangle::kUpper);
-      auto lower = (triangle == Triangle::kLower);
-      status = PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
-                                      n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                      n, n, c_ld, c_offset, c_buffer,
-                                      ConstantOne<T>(), program,
-                                      false, c_rotated, false, upper, lower, false);
-      if (ErrorIn(status)) { return status; }
-
-      // Successfully finished the computation
-      return StatusCode::kSuccess;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+  // Loads the program from the database
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+
+  // Determines whether or not temporary matrices are needed
+  auto a_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                   ab_rotated == false;
+  auto b_no_temp = ab_one == n_ceiled && ab_two == k_ceiled && b_ld == n_ceiled && b_offset == 0 &&
+                   ab_rotated == false;
+
+  // Creates the temporary matrices
+  auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto b_temp = (b_no_temp) ? b_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
+
+  // Events of all kernels (including pre/post processing kernels)
+  auto eventWaitList = std::vector<Event>();
+  auto emptyEventList = std::vector<Event>();
+
+  // Runs the pre-processing kernels. This transposes the matrices A and B, but also pads zeros to
+  // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
+  // case nothing has to be done, these kernels can be skipped.
+  if (!a_no_temp) {
+    auto eventProcessA = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
+                           ab_one, ab_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, false);
+    eventWaitList.push_back(eventProcessA);
+  }
+  if (!b_no_temp) {
+    auto eventProcessB = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessB.pointer(), emptyEventList,
+                           ab_one, ab_two, b_ld, b_offset, b_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, b_temp,
+                           ConstantOne<T>(), program,
+                           true, ab_rotated, false);
+    eventWaitList.push_back(eventProcessB);
+  }
+
+  // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
+  // modify the other triangle.
+  auto eventProcessC = Event();
+  PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
+                         n, n, c_ld, c_offset, c_buffer,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         ConstantOne<T>(), program,
+                         true, c_rotated, false);
+  eventWaitList.push_back(eventProcessC);
+
+  // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(n_ceiled));
+  kernel.SetArgument(1, static_cast<int>(k_ceiled));
+  kernel.SetArgument(2, GetRealArg(alpha));
+  kernel.SetArgument(3, GetRealArg(beta));
+  kernel.SetArgument(4, a_temp());
+  kernel.SetArgument(5, b_temp());
+  kernel.SetArgument(6, c_temp());
+
+  // Computes the global and local thread sizes
+  auto global = std::vector<size_t>{
+    (n_ceiled * db_["MDIMC"]) / db_["MWG"],
+    (n_ceiled * db_["NDIMC"]) / db_["NWG"]
+  };
+  auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
+
+  // Launches the kernel
+  auto eventKernel1 = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel1.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel1);
+
+  // Swaps the arguments for matrices A and B, and sets 'beta' to 1
+  auto one = static_cast<T>(1);
+  kernel.SetArgument(3, GetRealArg(one));
+  kernel.SetArgument(4, b_temp());
+  kernel.SetArgument(5, a_temp());
+
+  // Runs the kernel again
+  auto eventKernel2 = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel2.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel2);
+
+  // Runs the post-processing kernel
+  auto upper = (triangle == Triangle::kUpper);
+  auto lower = (triangle == Triangle::kLower);
+  PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         n, n, c_ld, c_offset, c_buffer,
+                         ConstantOne<T>(), program,
+                         false, c_rotated, false, upper, lower, false);
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xsyr2k.hpp b/src/routines/level3/xsyr2k.hpp
index 56185653..a02c6e16 100644
--- a/src/routines/level3/xsyr2k.hpp
+++ b/src/routines/level3/xsyr2k.hpp
@@ -30,13 +30,13 @@ class Xsyr2k: public Routine {
   Xsyr2k(Queue &queue, EventPointer event, const std::string &name = "SYR2K");
 
   // Templated-precision implementation of the routine
-  StatusCode DoSyr2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
-                     const size_t n, const size_t k,
-                     const T alpha,
-                     const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                     const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
-                     const T beta,
-                     const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoSyr2k(const Layout layout, const Triangle triangle, const Transpose ab_transpose,
+               const size_t n, const size_t k,
+               const T alpha,
+               const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+               const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld,
+               const T beta,
+               const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xsyrk.cpp b/src/routines/level3/xsyrk.cpp
index 438aa218..bd6c4b25 100644
--- a/src/routines/level3/xsyrk.cpp
+++ b/src/routines/level3/xsyrk.cpp
@@ -22,8 +22,7 @@ namespace clblast {
 // Constructor: forwards to base class constructor
 template <typename T>
 Xsyrk<T>::Xsyrk(Queue &queue, EventPointer event, const std::string &name):
-    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>()) {
-  source_string_ =
+    Routine(queue, event, name, {"Copy","Pad","Transpose","Padtranspose","Xgemm"}, PrecisionValue<T>(), {}, {
     #include "../../kernels/level3/level3.opencl"
     #include "../../kernels/level3/copy_fast.opencl"
     #include "../../kernels/level3/copy_pad.opencl"
@@ -32,14 +31,14 @@ Xsyrk<T>::Xsyrk(Queue &queue, EventPointer event, const std::string &name):
     #include "../../kernels/level3/xgemm_part1.opencl"
     #include "../../kernels/level3/xgemm_part2.opencl"
     #include "../../kernels/level3/xgemm_part3.opencl"
-  ;
+    }) {
 }
 
 // =================================================================================================
 
 // The main routine
 template <typename T>
-StatusCode Xsyrk<T>::DoSyrk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
+void Xsyrk<T>::DoSyrk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
                             const size_t n, const size_t k,
                             const T alpha,
                             const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
@@ -47,7 +46,7 @@ StatusCode Xsyrk<T>::DoSyrk(const Layout layout, const Triangle triangle, const
                             const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((n == 0) || (k == 0) ) { return StatusCode::kInvalidDimension; }
+  if ((n == 0) || (k == 0) ) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes whether or not the matrices are transposed in memory. This is based on their layout
   // (row or column-major) and whether or not they are requested to be pre-transposed.
@@ -65,10 +64,8 @@ StatusCode Xsyrk<T>::DoSyrk(const Layout layout, const Triangle triangle, const
   // space. Also tests that the leading dimensions of:
   //    matrix A cannot be less than N when rotated, or less than K when not-rotated
   //    matrix C cannot be less than N
-  auto status = TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
-  status = TestMatrixC(n, n, c_buffer, c_offset, c_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(a_one, a_two, a_buffer, a_offset, a_ld);
+  TestMatrixC(n, n, c_buffer, c_offset, c_ld);
 
   // Calculates the ceiled versions of n and k
   auto n_ceiled = Ceil(Ceil(n, db_["MWG"]), db_["NWG"]);
@@ -77,90 +74,76 @@ StatusCode Xsyrk<T>::DoSyrk(const Layout layout, const Triangle triangle, const
   // Decides which kernel to run: the upper-triangular or lower-triangular version
   auto kernel_name = (triangle == Triangle::kUpper) ? "XgemmUpper" : "XgemmLower";
 
-  // The padded/transposed input/output matrices: if memory allocation fails, throw an exception
-  try {
-
-    // Loads the program from the database
-    const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-
-    // Determines whether or not temporary matrices are needed
-    auto a_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
-                     a_rotated == false;
-
-    // Creates the temporary matrices
-    auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
-    auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
-
-    // Events of all kernels (including pre/post processing kernels)
-    auto eventWaitList = std::vector<Event>();
-    auto emptyEventList = std::vector<Event>();
-
-    // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
-    // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
-    // case nothing has to be done, these kernels can be skipped.
-    if (!a_no_temp) {
-      auto eventProcessA = Event();
-      status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
-                                      a_one, a_two, a_ld, a_offset, a_buffer,
-                                      n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
-                                      ConstantOne<T>(), program,
-                                      true, a_rotated, false);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventProcessA);
-    }
-
-    // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
-    // modify the other triangle.
-    auto eventProcessC = Event();
-    status = PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
-                                    n, n, c_ld, c_offset, c_buffer,
-                                    n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                    ConstantOne<T>(), program,
-                                    true, c_rotated, false);
-    if (ErrorIn(status)) { return status; }
-    eventWaitList.push_back(eventProcessC);
-
-    // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
-    try {
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the kernel arguments
-      kernel.SetArgument(0, static_cast<int>(n_ceiled));
-      kernel.SetArgument(1, static_cast<int>(k_ceiled));
-      kernel.SetArgument(2, GetRealArg(alpha));
-      kernel.SetArgument(3, GetRealArg(beta));
-      kernel.SetArgument(4, a_temp());
-      kernel.SetArgument(5, a_temp());
-      kernel.SetArgument(6, c_temp());
-
-      // Computes the global and local thread sizes
-      auto global = std::vector<size_t>{
-        (n_ceiled * db_["MDIMC"]) / db_["MWG"],
-        (n_ceiled * db_["NDIMC"]) / db_["NWG"]
-      };
-      auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
-
-      // Launches the kernel
-      auto eventKernel = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, eventKernel.pointer(), eventWaitList);
-      if (ErrorIn(status)) { return status; }
-      eventWaitList.push_back(eventKernel);
-
-      // Runs the post-processing kernel
-      auto upper = (triangle == Triangle::kUpper);
-      auto lower = (triangle == Triangle::kLower);
-      status = PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
-                                      n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
-                                      n, n, c_ld, c_offset, c_buffer,
-                                      ConstantOne<T>(), program,
-                                      false, c_rotated, false, upper, lower, false);
-      if (ErrorIn(status)) { return status; }
-
-
-      // Successfully finished the computation
-      return StatusCode::kSuccess;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+  // Loads the program from the database
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+
+  // Determines whether or not temporary matrices are needed
+  auto a_no_temp = a_one == n_ceiled && a_two == k_ceiled && a_ld == n_ceiled && a_offset == 0 &&
+                   a_rotated == false;
+
+  // Creates the temporary matrices
+  auto a_temp = (a_no_temp) ? a_buffer : Buffer<T>(context_, k_ceiled*n_ceiled);
+  auto c_temp = Buffer<T>(context_, n_ceiled*n_ceiled);
+
+  // Events of all kernels (including pre/post processing kernels)
+  auto eventWaitList = std::vector<Event>();
+  auto emptyEventList = std::vector<Event>();
+
+  // Runs the pre-processing kernel for matrix A. This transposes the matrix, but also pads zeros
+  // to fill it up until it reaches a certain multiple of size (kernel parameter dependent). In
+  // case nothing has to be done, these kernels can be skipped.
+  if (!a_no_temp) {
+    auto eventProcessA = Event();
+    PadCopyTransposeMatrix(queue_, device_, db_, eventProcessA.pointer(), emptyEventList,
+                           a_one, a_two, a_ld, a_offset, a_buffer,
+                           n_ceiled, k_ceiled, n_ceiled, 0, a_temp,
+                           ConstantOne<T>(), program,
+                           true, a_rotated, false);
+    eventWaitList.push_back(eventProcessA);
+  }
+
+  // Furthermore, also creates a (possibly padded) copy of matrix C, since it is not allowed to
+  // modify the other triangle.
+  auto eventProcessC = Event();
+  PadCopyTransposeMatrix(queue_, device_, db_, eventProcessC.pointer(), emptyEventList,
+                         n, n, c_ld, c_offset, c_buffer,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         ConstantOne<T>(), program,
+                         true, c_rotated, false);
+  eventWaitList.push_back(eventProcessC);
+
+  // Retrieves the XgemmUpper or XgemmLower kernel from the compiled binary
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the kernel arguments
+  kernel.SetArgument(0, static_cast<int>(n_ceiled));
+  kernel.SetArgument(1, static_cast<int>(k_ceiled));
+  kernel.SetArgument(2, GetRealArg(alpha));
+  kernel.SetArgument(3, GetRealArg(beta));
+  kernel.SetArgument(4, a_temp());
+  kernel.SetArgument(5, a_temp());
+  kernel.SetArgument(6, c_temp());
+
+  // Computes the global and local thread sizes
+  auto global = std::vector<size_t>{
+    (n_ceiled * db_["MDIMC"]) / db_["MWG"],
+    (n_ceiled * db_["NDIMC"]) / db_["NWG"]
+  };
+  auto local = std::vector<size_t>{db_["MDIMC"], db_["NDIMC"]};
+
+  // Launches the kernel
+  auto eventKernel = Event();
+  RunKernel(kernel, queue_, device_, global, local, eventKernel.pointer(), eventWaitList);
+  eventWaitList.push_back(eventKernel);
+
+  // Runs the post-processing kernel
+  auto upper = (triangle == Triangle::kUpper);
+  auto lower = (triangle == Triangle::kLower);
+  PadCopyTransposeMatrix(queue_, device_, db_, event_, eventWaitList,
+                         n_ceiled, n_ceiled, n_ceiled, 0, c_temp,
+                         n, n, c_ld, c_offset, c_buffer,
+                         ConstantOne<T>(), program,
+                         false, c_rotated, false, upper, lower, false);
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xsyrk.hpp b/src/routines/level3/xsyrk.hpp
index 7c075c26..de42b824 100644
--- a/src/routines/level3/xsyrk.hpp
+++ b/src/routines/level3/xsyrk.hpp
@@ -32,12 +32,12 @@ class Xsyrk: public Routine {
   Xsyrk(Queue &queue, EventPointer event, const std::string &name = "SYRK");
 
   // Templated-precision implementation of the routine
-  StatusCode DoSyrk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
-                    const size_t n, const size_t k,
-                    const T alpha,
-                    const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                    const T beta,
-                    const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
+  void DoSyrk(const Layout layout, const Triangle triangle, const Transpose a_transpose,
+              const size_t n, const size_t k,
+              const T alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const T beta,
+              const Buffer<T> &c_buffer, const size_t c_offset, const size_t c_ld);
 };
 
 // =================================================================================================
diff --git a/src/routines/level3/xtrmm.cpp b/src/routines/level3/xtrmm.cpp
index 74a82822..6bf77cfa 100644
--- a/src/routines/level3/xtrmm.cpp
+++ b/src/routines/level3/xtrmm.cpp
@@ -29,7 +29,7 @@ Xtrmm<T>::Xtrmm(Queue &queue, EventPointer event, const std::string &name):
 
 // The main routine
 template <typename T>
-StatusCode Xtrmm<T>::DoTrmm(const Layout layout, const Side side, const Triangle triangle,
+void Xtrmm<T>::DoTrmm(const Layout layout, const Side side, const Triangle triangle,
                             const Transpose a_transpose, const Diagonal diagonal,
                             const size_t m, const size_t n,
                             const T alpha,
@@ -37,15 +37,14 @@ StatusCode Xtrmm<T>::DoTrmm(const Layout layout, const Side side, const Triangle
                             const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld) {
 
   // Makes sure all dimensions are larger than zero
-  if ((m == 0) || (n == 0)) { return StatusCode::kInvalidDimension; }
+  if ((m == 0) || (n == 0)) { throw BLASError(StatusCode::kInvalidDimension); }
 
   // Computes the k dimension. This is based on whether or not matrix is A (on the left)
   // or B (on the right) in the Xgemm routine.
   auto k = (side == Side::kLeft) ? m : n;
 
   // Checks for validity of the triangular A matrix
-  auto status = TestMatrixA(k, k, a_buffer, a_offset, a_ld);
-  if (ErrorIn(status)) { return status; }
+  TestMatrixA(k, k, a_buffer, a_offset, a_ld);
 
   // Determines which kernel to run based on the layout (the Xgemm kernel assumes column-major as
   // default) and on whether we are dealing with an upper or lower triangle of the triangular matrix
@@ -57,74 +56,69 @@ StatusCode Xtrmm<T>::DoTrmm(const Layout layout, const Side side, const Triangle
   auto unit_diagonal = (diagonal == Diagonal::kUnit) ? true : false;
 
   // Temporary buffer for a copy of the triangular matrix
-  try {
-    auto temp_triangular = Buffer<T>(context_, k*k);
-
-    // Creates a general matrix from the triangular matrix to be able to run the regular Xgemm
-    // routine afterwards
+  auto temp_triangular = Buffer<T>(context_, k*k);
+
+  // Creates a general matrix from the triangular matrix to be able to run the regular Xgemm
+  // routine afterwards
+  const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
+  auto kernel = Kernel(program, kernel_name);
+
+  // Sets the arguments for the triangular-to-squared kernel
+  kernel.SetArgument(0, static_cast<int>(k));
+  kernel.SetArgument(1, static_cast<int>(a_ld));
+  kernel.SetArgument(2, static_cast<int>(a_offset));
+  kernel.SetArgument(3, a_buffer());
+  kernel.SetArgument(4, static_cast<int>(k));
+  kernel.SetArgument(5, static_cast<int>(k));
+  kernel.SetArgument(6, static_cast<int>(0));
+  kernel.SetArgument(7, temp_triangular());
+  kernel.SetArgument(8, static_cast<int>(unit_diagonal));
+
+  // Uses the common padding kernel's thread configuration. This is allowed, since the
+  // triangular-to-squared kernel uses the same parameters.
+  auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
+                                    Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
+  auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
+  auto kernelEvent = Event();
+  RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
+
+  // Synchronize now: 'DoGemm' does not accept a list of events to wait for
+  kernelEvent.WaitForCompletion();
+
+  // Runs the regular Xgemm code with either "B := alpha*A*B" or ...
+  if (side == Side::kLeft) {
+    DoGemm(layout, a_transpose, Transpose::kNo,
+           m, n, k,
+           alpha,
+           temp_triangular, 0, k,
+           b_buffer, b_offset, b_ld,
+           static_cast<T>(0.0),
+           b_buffer, b_offset, b_ld);
+  }
+
+  // ... with "B := alpha*B*A". Note that A and B are now reversed.
+  else {
     try {
-      const auto program = GetProgramFromCache(context_, PrecisionValue<T>(), routine_name_);
-      auto kernel = Kernel(program, kernel_name);
-
-      // Sets the arguments for the triangular-to-squared kernel
-      kernel.SetArgument(0, static_cast<int>(k));
-      kernel.SetArgument(1, static_cast<int>(a_ld));
-      kernel.SetArgument(2, static_cast<int>(a_offset));
-      kernel.SetArgument(3, a_buffer());
-      kernel.SetArgument(4, static_cast<int>(k));
-      kernel.SetArgument(5, static_cast<int>(k));
-      kernel.SetArgument(6, static_cast<int>(0));
-      kernel.SetArgument(7, temp_triangular());
-      kernel.SetArgument(8, static_cast<int>(unit_diagonal));
-
-      // Uses the common padding kernel's thread configuration. This is allowed, since the
-      // triangular-to-squared kernel uses the same parameters.
-      auto global = std::vector<size_t>{Ceil(CeilDiv(k, db_["PAD_WPTX"]), db_["PAD_DIMX"]),
-                                        Ceil(CeilDiv(k, db_["PAD_WPTY"]), db_["PAD_DIMY"])};
-      auto local = std::vector<size_t>{db_["PAD_DIMX"], db_["PAD_DIMY"]};
-      auto kernelEvent = Event();
-      status = RunKernel(kernel, queue_, device_, global, local, kernelEvent.pointer());
-      if (ErrorIn(status)) { return status; }
-
-      // Synchronize now: 'DoGemm' does not accept a list of events to wait for
-      kernelEvent.WaitForCompletion();
-
-      // Runs the regular Xgemm code with either "B := alpha*A*B" or ...
-      if (side == Side::kLeft) {
-        status = DoGemm(layout, a_transpose, Transpose::kNo,
-                        m, n, k,
-                        alpha,
-                        temp_triangular, 0, k,
-                        b_buffer, b_offset, b_ld,
-                        static_cast<T>(0.0),
-                        b_buffer, b_offset, b_ld);
-      }
-
-      // ... with "B := alpha*B*A". Note that A and B are now reversed.
-      else {
-        status = DoGemm(layout, Transpose::kNo, a_transpose,
-                        m, n, k,
-                        alpha,
-                        b_buffer, b_offset, b_ld,
-                        temp_triangular, 0, k,
-                        static_cast<T>(0.0),
-                        b_buffer, b_offset, b_ld);
-
-        // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
-        switch(status) {
-          case StatusCode::kInvalidMatrixA:      status = StatusCode::kInvalidMatrixB; break;
-          case StatusCode::kInvalidMatrixB:      status = StatusCode::kInvalidMatrixA; break;
-          case StatusCode::kInvalidLeadDimA:     status = StatusCode::kInvalidLeadDimB; break;
-          case StatusCode::kInvalidLeadDimB:     status = StatusCode::kInvalidLeadDimA; break;
-          case StatusCode::kInsufficientMemoryA: status = StatusCode::kInsufficientMemoryB; break;
-          case StatusCode::kInsufficientMemoryB: status = StatusCode::kInsufficientMemoryA; break;
-        }
+      DoGemm(layout, Transpose::kNo, a_transpose,
+             m, n, k,
+             alpha,
+             b_buffer, b_offset, b_ld,
+             temp_triangular, 0, k,
+             static_cast<T>(0.0),
+             b_buffer, b_offset, b_ld);
+    } catch (BLASError &e) {
+      // A and B are now reversed, so also reverse the error codes returned from the Xgemm routine
+      switch(e.status()) {
+        case StatusCode::kInvalidMatrixA:      throw BLASError(StatusCode::kInvalidMatrixB, e.details());
+        case StatusCode::kInvalidMatrixB:      throw BLASError(StatusCode::kInvalidMatrixA, e.details());
+        case StatusCode::kInvalidLeadDimA:     throw BLASError(StatusCode::kInvalidLeadDimB, e.details());
+        case StatusCode::kInvalidLeadDimB:     throw BLASError(StatusCode::kInvalidLeadDimA, e.details());
+        case StatusCode::kInsufficientMemoryA: throw BLASError(StatusCode::kInsufficientMemoryB, e.details());
+        case StatusCode::kInsufficientMemoryB: throw BLASError(StatusCode::kInsufficientMemoryA, e.details());
+        default:                               throw;
       }
-
-      // Return the status of the Xgemm routine
-      return status;
-    } catch (...) { return StatusCode::kInvalidKernel; }
-  } catch (...) { return StatusCode::kTempBufferAllocFailure; }
+    }
+  }
 }
 
 // =================================================================================================
diff --git a/src/routines/level3/xtrmm.hpp b/src/routines/level3/xtrmm.hpp
index 186a120e..967bf132 100644
--- a/src/routines/level3/xtrmm.hpp
+++ b/src/routines/level3/xtrmm.hpp
@@ -38,12 +38,12 @@ class Xtrmm: public Xgemm<T> {
   Xtrmm(Queue &queue, EventPointer event, const std::string &name = "TRMM");
 
   // Templated-precision implementation of the routine
-  StatusCode DoTrmm(const Layout layout, const Side side, const Triangle triangle,
-                    const Transpose a_transpose, const Diagonal diagonal,
-                    const size_t m, const size_t n,
-                    const T alpha,
-                    const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
-                    const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld);
+  void DoTrmm(const Layout layout, const Side side, const Triangle triangle,
+              const Transpose a_transpose, const Diagonal diagonal,
+              const size_t m, const size_t n,
+              const T alpha,
+              const Buffer<T> &a_buffer, const size_t a_offset, const size_t a_ld,
+              const Buffer<T> &b_buffer, const size_t b_offset, const size_t b_ld);
 };
 
 // =================================================================================================