Added data-preparation function for the TRSV tests and special nan/inf checks in the error checking to make the tests pass

4 files changed, 64 insertions, 2 deletions

diff --git a/src/utilities/utilities.cpp b/src/utilities/utilities.cpp
index fd198c0d..b2ed2f0c 100644
--- a/src/utilities/utilities.cpp
+++ b/src/utilities/utilities.cpp
@@ -18,6 +18,7 @@
 #include <chrono>
 #include <random>
 #include <iomanip>
+#include <cmath>
 
 namespace clblast {
 // =================================================================================================
@@ -118,6 +119,14 @@ double2 ConstantNegOne() {
   return {-1.0, 0.0};
 }
 
+// Returns the absolute value of a scalar
+template <typename T> T AbsoluteValue(const T value) { return std::fabs(value); }
+template float AbsoluteValue<float>(const float);
+template double AbsoluteValue<double>(const double);
+template <> half AbsoluteValue(const half value) { return FloatToHalf(std::fabs(HalfToFloat(value))); }
+template <> float2 AbsoluteValue(const float2 value) { return std::abs(value); }
+template <> double2 AbsoluteValue(const double2 value) { return std::abs(value); }
+
 // =================================================================================================
 
 // Implements the string conversion using std::to_string if possible
diff --git a/src/utilities/utilities.hpp b/src/utilities/utilities.hpp
index 757f1b5e..2c13658b 100644
--- a/src/utilities/utilities.hpp
+++ b/src/utilities/utilities.hpp
@@ -114,6 +114,10 @@ T ConstantOne();
 template <typename T>
 T ConstantNegOne();
 
+// Returns the absolute value of a scalar
+template <typename T>
+T AbsoluteValue(const T value);
+
 // =================================================================================================
 
 // Structure containing all possible arguments for test clients, including their default values
diff --git a/test/correctness/tester.cpp b/test/correctness/tester.cpp
index c0ed5ba4..dd2f3f99 100644
--- a/test/correctness/tester.cpp
+++ b/test/correctness/tester.cpp
@@ -40,6 +40,12 @@ float getAbsoluteErrorMargin<half>() {
   return 0.10f; // especially small values are inaccurate for half-precision
 }
 
+// Error margin: numbers beyond this value are considered equal to inf or NaN
+template <typename T>
+T getAlmostInfNumber() {
+  return T{1e35}; // used for correctness testing of TRSV and TRSM routines
+}
+
 // Maximum number of test results printed on a single line
 template <typename T, typename U> const size_t Tester<T,U>::kResultsPerLine = size_t{64};
 
@@ -426,8 +432,19 @@ bool TestSimilarityNear(const T val1, const T val2,
   if (val1 == val2) {
     return true;
   }
-  // Handles cases with both results NaN
-  else if (std::isnan(val1) && std::isnan(val2)) {
+  // Handles cases with both results NaN or inf
+  else if ((std::isnan(val1) && std::isnan(val2)) || (std::isinf(val1) && std::isinf(val2))) {
+    return true;
+  }
+  // Also considers it OK if one of the results in NaN and the other is inf
+  // Note: for TRSV and TRSM routines
+  else if ((std::isnan(val1) && std::isinf(val2)) || (std::isinf(val1) && std::isnan(val2))) {
+    return true;
+  }
+  // Also considers it OK if one of the values is super large and the other is inf or NaN
+  // Note: for TRSV and TRSM routines
+  else if ((std::abs(val1) > getAlmostInfNumber<T>() && (std::isinf(val2) || std::isnan(val2))) ||
+           (std::abs(val2) > getAlmostInfNumber<T>() && (std::isinf(val1) || std::isnan(val1)))) {
     return true;
   }
   // The values are zero or very small: the relative error is less meaningful
diff --git a/test/routines/level2/xtrsv.hpp b/test/routines/level2/xtrsv.hpp
index 67094b3d..811feac5 100644
--- a/test/routines/level2/xtrsv.hpp
+++ b/test/routines/level2/xtrsv.hpp
@@ -29,6 +29,35 @@
 namespace clblast {
 // =================================================================================================
 
+// Prepares the data
+template <typename T>
+void PrepareData(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+  if (args.a_ld < args.n) { return; }
+
+  // Copies input buffers to the host
+  std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
+  std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
+  buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);
+  buffers.x_vec.Read(queue, args.x_size, x_vec_cpu);
+
+  // Generates 'proper' input for the TRSV routine
+  // TODO: Improve this, currently loosely based on clBLAS's implementation
+  for (auto i = size_t{0}; i < args.n; ++i) {
+    auto diagonal = a_mat_cpu[i*args.a_ld + i + args.a_offset];
+    diagonal = AbsoluteValue(diagonal) + static_cast<T>(args.n / size_t{4});
+    for (auto j = size_t{0}; j < args.n; ++j) {
+      a_mat_cpu[j*args.a_ld + i + args.a_offset] /= T{2.0};
+    }
+    a_mat_cpu[i*args.a_ld + i + args.a_offset] = diagonal;
+    x_vec_cpu[i * args.x_inc + args.x_offset] /= T{2.0};
+  }
+
+  // Copies input buffers back to the OpenCL device
+  buffers.a_mat.Write(queue, args.a_size, a_mat_cpu);
+  buffers.x_vec.Write(queue, args.x_size, x_vec_cpu);
+  return;
+}
+
 // See comment at top of file for a description of the class
 template <typename T>
 class TestXtrsv {
@@ -71,6 +100,7 @@ class TestXtrsv {
 
   // Describes how to run the CLBlast routine
   static StatusCode RunRoutine(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+    PrepareData(args, buffers, queue);
     auto queue_plain = queue();
     auto event = cl_event{};
     auto status = Trsv<T>(args.layout, args.triangle, args.a_transpose, args.diagonal,
@@ -85,6 +115,7 @@ class TestXtrsv {
   // Describes how to run the clBLAS routine (for correctness/performance comparison)
   #ifdef CLBLAST_REF_CLBLAS
     static StatusCode RunReference1(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      PrepareData(args, buffers, queue);
       auto queue_plain = queue();
       auto event = cl_event{};
       auto status = clblasXtrsv<T>(convertToCLBLAS(args.layout),
@@ -103,6 +134,7 @@ class TestXtrsv {
   // Describes how to run the CPU BLAS routine (for correctness/performance comparison)
   #ifdef CLBLAST_REF_CBLAS
     static StatusCode RunReference2(const Arguments<T> &args, Buffers<T> &buffers, Queue &queue) {
+      PrepareData(args, buffers, queue);
       std::vector<T> a_mat_cpu(args.a_size, static_cast<T>(0));
       std::vector<T> x_vec_cpu(args.x_size, static_cast<T>(0));
       buffers.a_mat.Read(queue, args.a_size, a_mat_cpu);