DataHandlerGPU: Test core functionality (#21)

elsa/core/tests/CMakeLists.txt

@@ -12,7 +12,7 @@ ELSA_TEST(PartitionDescriptor)
 ELSA_TEST(RandomBlocksDescriptor)
 ELSA_TEST(DataContainer)
 ELSA_TEST(DataHandlers)
-ELSA_TEST(DataHandlerMapCPU)
+ELSA_TEST(DataHandlerMap)
 ELSA_TEST(LinearOperator)
 ELSA_TEST(ExpressionTemplates)
 

elsa/core/tests/test_BenchmarkExpressionTemplates.cpp

@@ -14,7 +14,11 @@
 #include <cstdlib>
 
 using namespace elsa;
-static const index_t dimension = 2;
+// use max of 256
+static const index_t dimension = 256;
+
+// dimension of memory critical benchmarks
+static const index_t dimensionMemCritic = 1024;
 
 TEST_CASE("Expression benchmark using Eigen with n=" + std::to_string(dimension) + "^3")
 {
@@ -85,18 +89,21 @@ TEST_CASE("Expression benchmark using expression templates with n=" + std::to_st
     {
         result = dc * dc2 - dc2 / dc3 + dc * dc3;
     };
+
+    BENCHMARK("reduction") { dc.sum(); };
 }
 
-TEST_CASE("Expression benchmark without expression templates with n=" + std::to_string(dimension)
+TEST_CASE("Expression benchmark using GPU expression templates with n=" + std::to_string(dimension)
           + "^3")
 {
     IndexVector_t numCoeff(3);
     numCoeff << dimension, dimension, dimension;
     DataDescriptor desc(numCoeff);
-    DataContainer dc(desc);
-    DataContainer dc2(desc);
-    DataContainer dc3(desc);
-    DataContainer result(desc);
+
+    DataContainer dc(desc, DataHandlerType::GPU);
+    DataContainer dc2(desc, DataHandlerType::GPU);
+    DataContainer dc3(desc, DataHandlerType::GPU);
+    DataContainer result(desc, DataHandlerType::GPU);
 
     for (index_t i = 0; i < dc.getSize(); ++i) {
         dc[i] = static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / 100.0));
@@ -104,9 +111,6 @@ TEST_CASE("Expression benchmark without expression templates with n=" + std::to_
         dc3[i] = static_cast<float>(rand()) / (static_cast<float>(RAND_MAX / 100.0));
     }
 
-    // to avoid using expression templates
-    using namespace elsa::detail;
-
     BENCHMARK("exp = dc - dc2;") { result = dc - dc2; };
 
     BENCHMARK("exp = dc - dc2 + dc;") { result = dc - dc2 + dc; };
@@ -119,4 +123,36 @@ TEST_CASE("Expression benchmark without expression templates with n=" + std::to_
     {
         result = dc * dc2 - dc2 / dc3 + dc * dc3;
     };
+
+    BENCHMARK("reduction") { dc.sum(); };
+}
+
+TEST_CASE("Expression benchmark using GPU expression templates with n="
+          + std::to_string(dimensionMemCritic) + "^3")
+{
+    index_t size = dimensionMemCritic * dimensionMemCritic * dimensionMemCritic;
+
+    Eigen::Matrix<float, Eigen::Dynamic, 1> randVec(size);
+    Eigen::Matrix<float, Eigen::Dynamic, 1> randVec2(size);
+
+    randVec.setRandom();
+    randVec2.setRandom();
+
+    IndexVector_t numCoeff(3);
+    numCoeff << dimensionMemCritic, dimensionMemCritic, dimensionMemCritic;
+    DataDescriptor desc(numCoeff);
+
+    DataContainer dc(desc, randVec, DataHandlerType::GPU);
+    DataContainer dc2(desc, randVec2, DataHandlerType::GPU);
+
+    BENCHMARK("GPU: dc2 = 1.2 * dc + dc2;") { dc2 = 1.2f * dc + dc2; };
+
+    BENCHMARK("Eigen: dc2 = 1.2 * dc + dc2;")
+    {
+        randVec2 = (randVec.array() * 1.2f + randVec2.array()).matrix();
+    };
+
+    BENCHMARK("GPU: reduction") { return dc.sum(); };
+
+    BENCHMARK("Eigen: reduction") { return randVec.sum(); };
 }

elsa/core/tests/test_DataHandlers.cpp

@@ -13,16 +13,35 @@
 #include "DataHandlerMapCPU.h"
 #include "testHelpers.h"
 
+#ifdef ELSA_CUDA_VECTOR
+#include "DataHandlerGPU.h"
+#include "DataHandlerMapGPU.h"
+#endif
+
 template <typename data_t>
 long elsa::useCount(const DataHandlerCPU<data_t>& dh)
 {
     return dh._data.use_count();
 }
 
+#ifdef ELSA_CUDA_VECTOR
+template <typename data_t>
+long elsa::useCount(const DataHandlerGPU<data_t>& dh)
+{
+    return dh._data.use_count();
+}
+#endif
+
 using namespace elsa;
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "", (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -92,9 +111,15 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Constructing DataHandler", "", (DataHandle
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing equality operator on DataHandler", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing equality operator on DataHandler", "",
                            (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -148,8 +173,14 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing equality operator on DataHandler",
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Assigning to DataHandlerCPU", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Assigning to DataHandlerCPU", "", (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -376,7 +407,12 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Assigning to DataHandlerCPU", "", (DataHan
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_TEST_CASE("Scenario: Cloning DataHandler", "", DataHandlerCPU<float>,
+                   DataHandlerGPU<float>)
+#else
 TEMPLATE_TEST_CASE("Scenario: Cloning DataHandler", "", DataHandlerCPU<float>)
+#endif
 {
     GIVEN("some DataHandler")
     {
@@ -403,9 +439,15 @@ TEMPLATE_TEST_CASE("Scenario: Cloning DataHandler", "", DataHandlerCPU<float>)
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the reduction operatios of DataHandler", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the reduction operatios of DataHandler", "",
                            (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -421,7 +463,7 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the reduction operatios of DataHan
 
             THEN("the reductions work as expected")
             {
-                REQUIRE(dh.sum() == randVec.sum());
+                REQUIRE(checkSameNumbers(dh.sum(), randVec.sum()));
                 REQUIRE(dh.l1Norm() == Approx(randVec.array().abs().sum()));
                 REQUIRE(dh.lInfNorm() == Approx(randVec.array().abs().maxCoeff()));
                 REQUIRE(dh.squaredL2Norm() == Approx(randVec.squaredNorm()));
@@ -450,9 +492,15 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the reduction operatios of DataHan
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the element-wise operations of DataHandler", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the element-wise operations of DataHandler", "",
                            (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -504,24 +552,26 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the element-wise operations of Dat
                 dh = oldDh;
                 dh *= dh2;
                 for (index_t i = 0; i < size; ++i)
-                    REQUIRE(dh[i] == oldDh[i] * dh2[i]);
+                    REQUIRE(checkSameNumbers(dh[i], oldDh[i] * dh2[i]));
 
                 dh = oldDh;
                 dh *= *dhMap;
                 for (index_t i = 0; i < size; ++i)
-                    REQUIRE(dh[i] == oldDh[i] * dh2[i]);
+                    REQUIRE(checkSameNumbers(dh[i], oldDh[i] * dh2[i]));
 
                 dh = oldDh;
                 dh /= dh2;
                 for (index_t i = 0; i < size; ++i)
                     if (dh2[i] != data_t(0))
-                        REQUIRE(checkSameNumbers(dh[i], oldDh[i] / dh2[i]));
+                        // due to floating point arithmetic less precision
+                        REQUIRE(checkSameNumbers(dh[i], oldDh[i] / dh2[i], 100));
 
                 dh = oldDh;
                 dh /= *dhMap;
                 for (index_t i = 0; i < size; ++i)
                     if (dh2[i] != data_t(0))
-                        REQUIRE(checkSameNumbers(dh[i], oldDh[i] / dh2[i]));
+                        // due to floating point arithmetic less precision
+                        REQUIRE(checkSameNumbers(dh[i], oldDh[i] / dh2[i], 100));
             }
 
             THEN("the element-wise binary scalar operations work as expected")
@@ -561,8 +611,14 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the element-wise operations of Dat
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Referencing blocks of DataHandler", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Referencing blocks of DataHandler", "", (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 
@@ -601,8 +657,14 @@ TEMPLATE_PRODUCT_TEST_CASE("Scenario: Referencing blocks of DataHandler", "", (D
     }
 }
 
+#ifdef ELSA_CUDA_VECTOR
+TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the copy-on-write mechanism", "",
+                           (DataHandlerCPU, DataHandlerGPU),
+                           (float, double, std::complex<float>, std::complex<double>, index_t))
+#else
 TEMPLATE_PRODUCT_TEST_CASE("Scenario: Testing the copy-on-write mechanism", "", (DataHandlerCPU),
                            (float, double, std::complex<float>, std::complex<double>, index_t))
+#endif
 {
     using data_t = typename TestType::value_type;
 

elsa/core/tests/test_elsaDefines.cpp

@@ -7,6 +7,7 @@
  */
 
 #include <catch2/catch.hpp>
+#include <iostream>
 #include "elsaDefines.h"
 
 using namespace elsa;
@@ -36,4 +37,16 @@ SCENARIO("Testing compile-time predicates")
     static_assert(!std::is_same_v<float, GetFloatingPointType_t<double>>);
 
     REQUIRE(true);
+}
+
+SCENARIO("Printing default handler type")
+{
+    switch (defaultHandlerType) {
+        case DataHandlerType::CPU:
+            std::cout << "CPU" << std::endl;
+            break;
+        case DataHandlerType::GPU:
+            std::cout << "GPU" << std::endl;
+            break;
+    }
 }
\ No newline at end of file

elsa/test_routines/testHelpers.h

 #pragma once
 
 #include <type_traits>
-#include <complex.h>
+#include <complex>
+#include "elsaDefines.h"
 
 /**
  * \brief comparing two number types for approximate equality for complex and regular number
@@ -13,15 +14,22 @@
  * The CHECK(...) assertion in the function ensures that the values are reported when the test fails
  */
 template <typename T>
-bool checkSameNumbers(T right, T left)
+bool checkSameNumbers(T left, T right, int epsilonFactor = 1)
 {
+    using numericalBaseType = elsa::GetFloatingPointType_t<T>;
+
+    numericalBaseType eps = std::numeric_limits<numericalBaseType>::epsilon()
+                            * static_cast<numericalBaseType>(epsilonFactor)
+                            * static_cast<numericalBaseType>(100);
+
     if constexpr (std::is_same_v<T,
                                  std::complex<float>> || std::is_same_v<T, std::complex<double>>) {
-        CHECK(Approx(right.real()) == left.real());
-        CHECK(Approx(right.imag()) == left.imag());
-        return Approx(right.real()) == left.real() && Approx(right.imag()) == left.imag();
+        CHECK(Approx(left.real()).epsilon(eps) == right.real());
+        CHECK(Approx(left.imag()).epsilon(eps) == right.imag());
+        return Approx(left.real()).epsilon(eps) == right.real()
+               && Approx(left.imag()).epsilon(eps) == right.imag();
     } else {
-        CHECK(Approx(right) == left);
-        return Approx(right) == left;
+        CHECK(Approx(left).epsilon(eps) == right);
+        return Approx(left).epsilon(eps) == right;
     }
 }