Added initial version of GEMV including tester and performance client

6 files changed, 368 insertions, 17 deletions

diff --git a/test/correctness/routines/xgemv.cc b/test/correctness/routines/xgemv.cc
new file mode 100644
index 00000000..94ae147f
--- /dev/null
+++ b/test/correctness/routines/xgemv.cc
@@ -0,0 +1,94 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under the MIT license. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file implements the tests for the Xgemv routine. It is based on the TestAXY class.
+//
+// =================================================================================================
+
+#include "wrapper_clblas.h"
+#include "correctness/testaxy.h"
+
+namespace clblast {
+// =================================================================================================
+
+// The correctness tester, containing the function calls to CLBlast and to clBLAS for comparison.
+template <typename T>
+void XgemvTest(int argc, char *argv[], const bool silent, const std::string &name) {
+
+  // Creates the CLBlast lambda
+  auto clblast_lambda = [](const Arguments<T> &args,
+                           const Buffer &a_mat, const Buffer &x_vec, const Buffer &y_vec,
+                           CommandQueue &queue) -> StatusCode {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    return Gemv(args.layout, args.a_transpose, args.m, args.n, args.alpha,
+                a_mat(), args.a_offset, args.a_ld,
+                x_vec(), args.x_offset, args.x_inc, args.beta,
+                y_vec(), args.y_offset, args.y_inc,
+                &queue_plain, &event);
+  };
+
+  // Creates the clBLAS lambda (for comparison)
+  auto clblas_lambda = [](const Arguments<T> &args,
+                          const Buffer &a_mat, const Buffer &x_vec, const Buffer &y_vec,
+                          CommandQueue &queue) -> StatusCode {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXgemv(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasTranspose>(args.a_transpose),
+                              args.m, args.n, args.alpha,
+                              a_mat(), args.a_offset, args.a_ld,
+                              x_vec(), args.x_offset, args.x_inc, args.beta,
+                              y_vec(), args.y_offset, args.y_inc,
+                              1, &queue_plain, 0, nullptr, &event);
+    return static_cast<StatusCode>(status);
+  };
+
+  // Selects the platform and device on which to test (command-line options)
+  auto help = std::string{"Options given/available:\n"};
+  const auto platform_id = GetArgument(argc, argv, help, kArgPlatform, size_t{0});
+  const auto device_id = GetArgument(argc, argv, help, kArgDevice, size_t{0});
+  if (!silent) { fprintf(stdout, "\n* %s\n", help.c_str()); }
+
+  // Initializes the other arguments relevant for this routine
+  auto args = Arguments<T>{};
+  const auto options = std::vector<std::string>{kArgM, kArgN, kArgLayout, kArgATransp,
+                                                kArgALeadDim, kArgXInc, kArgYInc,
+                                                kArgAOffset, kArgXOffset, kArgYOffset};
+
+  // Creates a tester
+  TestAXY<T> tester{platform_id, device_id, name, options, clblast_lambda, clblas_lambda};
+
+  // Loops over the test-cases from a data-layout point of view
+  for (auto &layout: {Layout::kRowMajor, Layout::kColMajor}) {
+    args.layout = layout;
+    for (auto &a_transpose: {Transpose::kNo, Transpose::kYes}) {
+      args.a_transpose = a_transpose;
+      const auto case_name = ToString(layout)+" "+ToString(a_transpose);
+
+      // Runs the tests
+      tester.TestRegular(args, case_name);
+      tester.TestInvalidBufferSizes(args, case_name);
+    }
+  }
+}
+
+// =================================================================================================
+} // namespace clblast
+
+// Main function (not within the clblast namespace)
+int main(int argc, char *argv[]) {
+  clblast::XgemvTest<float>(argc, argv, false, "SGEMV");
+  //clblast::XgemvTest<double>(argc, argv, true, "DGEMV");
+  //clblast::XgemvTest<clblast::float2>(argc, argv, true, "CGEMV");
+  //clblast::XgemvTest<clblast::double2>(argc, argv, true, "ZGEMV");
+  return 0;
+}
+
+// =================================================================================================
diff --git a/test/correctness/testaxy.cc b/test/correctness/testaxy.cc
index 1e01a0e8..ed0b06ab 100644
--- a/test/correctness/testaxy.cc
+++ b/test/correctness/testaxy.cc
@@ -49,12 +49,6 @@ template <typename T>
 void TestAXY<T>::TestRegular(Arguments<T> &args, const std::string &name) {
   TestStart("regular behaviour", name);
 
-  // Computes whether or not the matrix is transposed. Note that we assume a default of
-  // column-major and no-transpose. If one of them is different (but not both), then rotated
-  // is considered true.
-  auto a_rotated = (args.layout == Layout::kColMajor && args.a_transpose != Transpose::kNo) ||
-                   (args.layout == Layout::kRowMajor && args.a_transpose == Transpose::kNo);
-
   // Iterates over the dimension for the matrix and vectors
   for (auto &m: kMatrixVectorDims) {
     args.m = m;
@@ -62,7 +56,12 @@ void TestAXY<T>::TestRegular(Arguments<T> &args, const std::string &name) {
       args.n = n;
 
       // Computes the second dimension of the matrix taking the rotation into account
-      auto a_two = (a_rotated) ? m : n;
+      auto a_two = (args.layout == Layout::kRowMajor) ? args.m : args.n;
+
+      // Computes the vector sizes in case the matrix is transposed
+      auto a_transposed = (args.a_transpose == Transpose::kYes);
+      auto m_real = (a_transposed) ? n : m;
+      auto n_real = (a_transposed) ? m : n;
 
       // Iterates over the leading-dimension values and the offsets of the matrix
       for (auto &a_ld: kMatrixVectorDims) {
@@ -82,8 +81,8 @@ void TestAXY<T>::TestRegular(Arguments<T> &args, const std::string &name) {
 
                   // Computes the buffer sizes
                   auto a_size = a_two * a_ld + a_offset;
-                  auto x_size = n * x_inc + x_offset;
-                  auto y_size = n * y_inc + y_offset;
+                  auto x_size = n_real * x_inc + x_offset;
+                  auto y_size = m_real * y_inc + y_offset;
                   if (a_size < 1 || x_size < 1 || y_size < 1) { continue; }
 
                   // Creates the OpenCL buffers
@@ -124,15 +123,15 @@ void TestAXY<T>::TestRegular(Arguments<T> &args, const std::string &name) {
 
                       // Checks for differences in the output
                       auto errors = size_t{0};
-                      for (auto idn=size_t{0}; idn<n; ++idn) {
-                        auto index = idn*y_inc + y_offset;
+                      for (auto idm=size_t{0}; idm<m_real; ++idm) {
+                        auto index = idm*y_inc + y_offset;
                         if (!TestSimilarity(r_result[index], s_result[index], kErrorMargin)) {
                           errors++;
                         }
                       }
 
                       // Tests the error count (should be zero)
-                      TestErrorCount(errors, n, args);
+                      TestErrorCount(errors, m_real, args);
                     }
                   }
                 }
@@ -158,6 +157,9 @@ void TestAXY<T>::TestInvalidBufferSizes(Arguments<T> &args, const std::string &n
   args.m = kBufferSize;
   args.n = kBufferSize;
   args.a_ld = kBufferSize;
+  args.a_offset = 0;
+  args.x_offset = 0;
+  args.y_offset = 0;
 
   // Iterates over test buffer sizes
   const std::vector<size_t> kMatrixSizes = {0, kBufferSize*kBufferSize-1, kBufferSize*kBufferSize};
diff --git a/test/performance/client.cc b/test/performance/client.cc
index ddaea0e1..3b07970c 100644
--- a/test/performance/client.cc
+++ b/test/performance/client.cc
@@ -26,8 +26,12 @@ template <typename T>
 void ClientXY(int argc, char *argv[], Routine2<T> client_routine,
               const std::vector<std::string> &options) {
 
+  // Function to determine how to find the default value of the leading dimension of matrix A.
+  // Note: this is not relevant for this client but given anyway.
+  auto default_ld_a = [](const Arguments<T> args) { return args.n; };
+
   // Simple command line argument parser with defaults
-  auto args = ParseArguments<T>(argc, argv, options);
+  auto args = ParseArguments<T>(argc, argv, options, default_ld_a);
   if (args.print_help) { return; }
 
   // Prints the header of the output table
@@ -81,13 +85,94 @@ template void ClientXY<double2>(int, char **, Routine2<double2>, const std::vect
 
 // =================================================================================================
 
+// This is the matrix-vector-vector variant of the set-up/tear-down client routine.
+template <typename T>
+void ClientAXY(int argc, char *argv[], Routine3<T> client_routine,
+               const std::vector<std::string> &options) {
+
+  // Function to determine how to find the default value of the leading dimension of matrix A
+  auto default_ld_a = [](const Arguments<T> args) { return args.n; };
+
+  // Simple command line argument parser with defaults
+  auto args = ParseArguments<T>(argc, argv, options, default_ld_a);
+  if (args.print_help) { return; }
+
+  // Prints the header of the output table
+  PrintTableHeader(args.silent, options);
+
+  // Initializes OpenCL and the libraries
+  auto platform = Platform(args.platform_id);
+  auto device = Device(platform, kDeviceType, args.device_id);
+  auto context = Context(device);
+  auto queue = CommandQueue(context, device);
+  if (args.compare_clblas) { clblasSetup(); }
+
+  // Iterates over all "num_step" values jumping by "step" each time
+  auto s = size_t{0};
+  while(true) {
+
+    // Computes the second dimension of the matrix taking the rotation into account
+    auto a_two = (args.layout == Layout::kRowMajor) ? args.m : args.n;
+
+    // Computes the vector sizes in case the matrix is transposed
+    auto a_transposed = (args.a_transpose == Transpose::kYes);
+    auto m_real = (a_transposed) ? args.n : args.m;
+    auto n_real = (a_transposed) ? args.m : args.n;
+
+    // Computes the data sizes
+    auto a_size = a_two * args.a_ld + args.a_offset;
+    auto x_size = n_real*args.x_inc + args.x_offset;
+    auto y_size = m_real*args.y_inc + args.y_offset;
+
+    // Populates input host vectors with random data
+    std::vector<T> a_source(a_size);
+    std::vector<T> x_source(x_size);
+    std::vector<T> y_source(y_size);
+    PopulateVector(a_source);
+    PopulateVector(x_source);
+    PopulateVector(y_source);
+
+    // Creates the vectors on the device
+    auto a_buffer = Buffer(context, CL_MEM_READ_WRITE, a_size*sizeof(T));
+    auto x_buffer = Buffer(context, CL_MEM_READ_WRITE, x_size*sizeof(T));
+    auto y_buffer = Buffer(context, CL_MEM_READ_WRITE, y_size*sizeof(T));
+    a_buffer.WriteBuffer(queue, a_size*sizeof(T), a_source);
+    x_buffer.WriteBuffer(queue, x_size*sizeof(T), x_source);
+    y_buffer.WriteBuffer(queue, y_size*sizeof(T), y_source);
+
+    // Runs the routine-specific code
+    client_routine(args, a_buffer, x_buffer, y_buffer, queue);
+
+    // Makes the jump to the next step
+    ++s;
+    if (s >= args.num_steps) { break; }
+    args.m += args.step;
+    args.n += args.step;
+    args.a_ld += args.step;
+  }
+
+  // Cleans-up and returns
+  if (args.compare_clblas) { clblasTeardown(); }
+}
+
+// Compiles the above function
+template void ClientAXY<float>(int, char **, Routine3<float>, const std::vector<std::string>&);
+template void ClientAXY<double>(int, char **, Routine3<double>, const std::vector<std::string>&);
+template void ClientAXY<float2>(int, char **, Routine3<float2>, const std::vector<std::string>&);
+template void ClientAXY<double2>(int, char **, Routine3<double2>, const std::vector<std::string>&);
+
+// =================================================================================================
+
 // This is the matrix-matrix-matrix variant of the set-up/tear-down client routine.
 template <typename T>
 void ClientABC(int argc, char *argv[], Routine3<T> client_routine,
                      const std::vector<std::string> &options) {
 
+  // Function to determine how to find the default value of the leading dimension of matrix A
+  auto default_ld_a = [](const Arguments<T> args) { return args.m; };
+
   // Simple command line argument parser with defaults
-  auto args = ParseArguments<T>(argc, argv, options);
+  auto args = ParseArguments<T>(argc, argv, options, default_ld_a);
   if (args.print_help) { return; }
 
   // Prints the header of the output table
@@ -167,7 +252,8 @@ template void ClientABC<double2>(int, char **, Routine3<double2>, const std::vec
 // applicable, but are searched for anyway to be able to create one common argument parser. All
 // arguments have a default value in case they are not found.
 template <typename T>
-Arguments<T> ParseArguments(int argc, char *argv[], const std::vector<std::string> &options) {
+Arguments<T> ParseArguments(int argc, char *argv[], const std::vector<std::string> &options,
+                            const std::function<size_t(const Arguments<T>)> default_ld_a) {
   auto args = Arguments<T>{};
   auto help = std::string{"Options given/available:\n"};
 
@@ -193,7 +279,7 @@ Arguments<T> ParseArguments(int argc, char *argv[], const std::vector<std::strin
     if (o == kArgYOffset) { args.y_offset = GetArgument(argc, argv, help, kArgYOffset, size_t{0}); }
 
     // Matrix arguments
-    if (o == kArgALeadDim) { args.a_ld     = GetArgument(argc, argv, help, kArgALeadDim, args.k); }
+    if (o == kArgALeadDim) { args.a_ld     = GetArgument(argc, argv, help, kArgALeadDim, default_ld_a(args)); }
     if (o == kArgBLeadDim) { args.b_ld     = GetArgument(argc, argv, help, kArgBLeadDim, args.n); }
     if (o == kArgCLeadDim) { args.c_ld     = GetArgument(argc, argv, help, kArgCLeadDim, args.n); }
     if (o == kArgAOffset)  { args.a_offset = GetArgument(argc, argv, help, kArgAOffset, size_t{0}); }
diff --git a/test/performance/client.h b/test/performance/client.h
index 2b9991fe..5125844a 100644
--- a/test/performance/client.h
+++ b/test/performance/client.h
@@ -49,6 +49,9 @@ template <typename T>
 void ClientXY(int argc, char *argv[], Routine2<T> client_routine,
               const std::vector<std::string> &options);
 template <typename T>
+void ClientAXY(int argc, char *argv[], Routine3<T> client_routine,
+               const std::vector<std::string> &options);
+template <typename T>
 void ClientABC(int argc, char *argv[], Routine3<T> client_routine,
                const std::vector<std::string> &options);
 
@@ -57,7 +60,8 @@ void ClientABC(int argc, char *argv[], Routine3<T> client_routine,
 // Parses all command-line arguments, filling in the arguments structure. If no command-line
 // argument is given for a particular argument, it is filled in with a default value.
 template <typename T>
-Arguments<T> ParseArguments(int argc, char *argv[], const std::vector<std::string> &options);
+Arguments<T> ParseArguments(int argc, char *argv[], const std::vector<std::string> &options,
+                            const std::function<size_t(const Arguments<T>)> default_ld_a);
 
 // Retrieves only the precision command-line argument, since the above function is templated based
 // on the precision
diff --git a/test/performance/routines/xgemv.cc b/test/performance/routines/xgemv.cc
new file mode 100644
index 00000000..43222396
--- /dev/null
+++ b/test/performance/routines/xgemv.cc
@@ -0,0 +1,107 @@
+
+// =================================================================================================
+// This file is part of the CLBlast project. The project is licensed under Apache Version 2.0. This
+// project loosely follows the Google C++ styleguide and uses a tab-size of two spaces and a max-
+// width of 100 characters per line.
+//
+// Author(s):
+//   Cedric Nugteren <www.cedricnugteren.nl>
+//
+// This file implements the Xgemv command-line interface tester.
+//
+// =================================================================================================
+
+#include <string>
+#include <vector>
+#include <exception>
+
+#include "wrapper_clblas.h"
+#include "performance/client.h"
+
+namespace clblast {
+// =================================================================================================
+
+// The client, used for performance testing. It contains the function calls to CLBlast and to other
+// libraries to compare against.
+template <typename T>
+void PerformanceXgemv(const Arguments<T> &args,
+                      const Buffer &a_mat, const Buffer &x_vec, const Buffer &y_vec,
+                      CommandQueue &queue) {
+
+  // Creates the CLBlast lambda
+  auto clblast_lambda = [&args, &a_mat, &x_vec, &y_vec, &queue]() {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = Gemv(args.layout, args.a_transpose, args.m, args.n, args.alpha,
+                       a_mat(), args.a_offset, args.a_ld,
+                       x_vec(), args.x_offset, args.x_inc, args.beta,
+                       y_vec(), args.y_offset, args.y_inc,
+                       &queue_plain, &event);
+    clWaitForEvents(1, &event);
+    if (status != StatusCode::kSuccess) {
+      throw std::runtime_error("CLBlast error: "+ToString(static_cast<int>(status)));
+    }
+  };
+
+  // Creates the clBLAS lambda (for comparison)
+  auto clblas_lambda = [&args, &a_mat, &x_vec, &y_vec, &queue]() {
+    auto queue_plain = queue();
+    auto event = cl_event{};
+    auto status = clblasXgemv(static_cast<clblasOrder>(args.layout),
+                              static_cast<clblasTranspose>(args.a_transpose),
+                              args.m, args.n, args.alpha,
+                              a_mat(), args.a_offset, args.a_ld,
+                              x_vec(), args.x_offset, args.x_inc, args.beta,
+                              y_vec(), args.y_offset, args.y_inc,
+                              1, &queue_plain, 0, nullptr, &event);
+    clWaitForEvents(1, &event);
+    if (status != CL_SUCCESS) {
+      throw std::runtime_error("clBLAS error: "+ToString(static_cast<int>(status)));
+    }
+  };
+
+  // Runs the routines and collect the timings
+  auto ms_clblast = TimedExecution(args.num_runs, clblast_lambda);
+  auto ms_clblas = TimedExecution(args.num_runs, clblas_lambda);
+
+  // Prints the performance of both libraries
+  const auto flops = 2 * args.m * args.n;
+  const auto bytes = (args.m*args.n + 2*args.m + args.n) * sizeof(T);
+  const auto output_ints = std::vector<size_t>{args.m, args.n,
+                                               static_cast<size_t>(args.layout),
+                                               static_cast<size_t>(args.a_transpose),
+                                               args.a_ld, args.x_inc, args.y_inc,
+                                               args.a_offset, args.x_offset, args.y_offset};
+  const auto output_strings = std::vector<std::string>{ToString(args.alpha),
+                                                       ToString(args.beta)};
+  PrintTableRow(output_ints, output_strings, args.no_abbrv,
+                ms_clblast, ms_clblas, flops, bytes);
+}
+
+// =================================================================================================
+
+// Main function which calls the common client code with the routine-specific function as argument.
+void ClientXgemv(int argc, char *argv[]) {
+  const auto o = std::vector<std::string>{kArgM, kArgN, kArgLayout, kArgATransp,
+                                          kArgALeadDim, kArgXInc, kArgYInc,
+                                          kArgAOffset, kArgXOffset, kArgYOffset,
+                                          kArgAlpha, kArgBeta};
+  switch(GetPrecision(argc, argv)) {
+    case Precision::kHalf: throw std::runtime_error("Unsupported precision mode");
+    case Precision::kSingle: ClientAXY<float>(argc, argv, PerformanceXgemv<float>, o); break;
+    case Precision::kDouble: ClientAXY<double>(argc, argv, PerformanceXgemv<double>, o); break;
+    case Precision::kComplexSingle: ClientAXY<float2>(argc, argv, PerformanceXgemv<float2>, o); break;
+    case Precision::kComplexDouble: ClientAXY<double2>(argc, argv, PerformanceXgemv<double2>, o); break;
+  }
+}
+
+// =================================================================================================
+} // namespace clblast
+
+// Main function (not within the clblast namespace)
+int main(int argc, char *argv[]) {
+  clblast::ClientXgemv(argc, argv);
+  return 0;
+}
+
+// =================================================================================================
diff --git a/test/wrapper_clblas.h b/test/wrapper_clblas.h
index 7c71fcaa..093a8742 100644
--- a/test/wrapper_clblas.h
+++ b/test/wrapper_clblas.h
@@ -74,6 +74,64 @@ clblasStatus clblasXaxpy(
 // =================================================================================================
 // BLAS level-2 (matrix-vector) routines
 
+// Calls {clblasSgemv, clblasDgemv, clblasCgemv, clblasZgemv} with the arguments forwarded.
+clblasStatus clblasXgemv(
+  clblasOrder layout, clblasTranspose tran_a, size_t m, size_t n, float alpha,
+  const cl_mem a_mat, size_t a_offset, size_t a_ld,
+  const cl_mem x_vec, size_t x_offset, size_t x_inc, float beta,
+  const cl_mem y_vec, size_t y_offset, size_t y_inc,
+  cl_uint num_queues, cl_command_queue *queues,
+  cl_uint num_wait_events, const cl_event *wait_events, cl_event *events) {
+    return clblasSgemv(layout, tran_a, m, n, alpha,
+                       a_mat, a_offset, a_ld,
+                       x_vec, x_offset, static_cast<int>(x_inc), beta,
+                       y_vec, y_offset, static_cast<int>(y_inc),
+                       num_queues, queues, num_wait_events, wait_events, events);
+}
+clblasStatus clblasXgemv(
+  clblasOrder layout, clblasTranspose tran_a, size_t m, size_t n, double alpha,
+  const cl_mem a_mat, size_t a_offset, size_t a_ld,
+  const cl_mem x_vec, size_t x_offset, size_t x_inc, double beta,
+  const cl_mem y_vec, size_t y_offset, size_t y_inc,
+  cl_uint num_queues, cl_command_queue *queues,
+  cl_uint num_wait_events, const cl_event *wait_events, cl_event *events) {
+    return clblasDgemv(layout, tran_a, m, n, alpha,
+                       a_mat, a_offset, a_ld,
+                       x_vec, x_offset, static_cast<int>(x_inc), beta,
+                       y_vec, y_offset, static_cast<int>(y_inc),
+                       num_queues, queues, num_wait_events, wait_events, events);
+}
+clblasStatus clblasXgemv(
+  clblasOrder layout, clblasTranspose tran_a, size_t m, size_t n, float2 alpha,
+  const cl_mem a_mat, size_t a_offset, size_t a_ld,
+  const cl_mem x_vec, size_t x_offset, size_t x_inc, float2 beta,
+  const cl_mem y_vec, size_t y_offset, size_t y_inc,
+  cl_uint num_queues, cl_command_queue *queues,
+  cl_uint num_wait_events, const cl_event *wait_events, cl_event *events) {
+    auto cl_alpha = cl_float2{{alpha.real(), alpha.imag()}};
+    auto cl_beta = cl_float2{{beta.real(), beta.imag()}};
+    return clblasCgemv(layout, tran_a, m, n, cl_alpha,
+                       a_mat, a_offset, a_ld,
+                       x_vec, x_offset, static_cast<int>(x_inc), cl_beta,
+                       y_vec, y_offset, static_cast<int>(y_inc),
+                       num_queues, queues, num_wait_events, wait_events, events);
+}
+clblasStatus clblasXgemv(
+  clblasOrder layout, clblasTranspose tran_a, size_t m, size_t n, double2 alpha,
+  const cl_mem a_mat, size_t a_offset, size_t a_ld,
+  const cl_mem x_vec, size_t x_offset, size_t x_inc, double2 beta,
+  const cl_mem y_vec, size_t y_offset, size_t y_inc,
+  cl_uint num_queues, cl_command_queue *queues,
+  cl_uint num_wait_events, const cl_event *wait_events, cl_event *events) {
+    auto cl_alpha = cl_double2{{alpha.real(), alpha.imag()}};
+    auto cl_beta = cl_double2{{beta.real(), beta.imag()}};
+    return clblasZgemv(layout, tran_a, m, n, cl_alpha,
+                       a_mat, a_offset, a_ld,
+                       x_vec, x_offset, static_cast<int>(x_inc), cl_beta,
+                       y_vec, y_offset, static_cast<int>(y_inc),
+                       num_queues, queues, num_wait_events, wait_events, events);
+}
+
 // =================================================================================================
 // BLAS level-3 (matrix-matrix) routines