Fixed isfinite implementation and added test

core/test/base/CMakeLists.txt

@@ -9,6 +9,7 @@ ginkgo_create_test(extended_float)
 ginkgo_create_test(executor)
 ginkgo_create_test(iterator_factory)
 ginkgo_create_test(lin_op)
+ginkgo_create_test(math)
 ginkgo_create_test(matrix_data)
 ginkgo_create_test(mtx_io)
 ginkgo_create_test(perturbation)

core/test/base/math.cpp

@@ -0,0 +1,107 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2019, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/base/math.hpp>
+
+
+#include <cmath>
+#include <complex>
+#include <limits>
+
+
+#include <gtest/gtest.h>
+
+
+namespace {
+
+
+template <typename T>
+void test_real_isfinite()
+{
+    using limits = std::numeric_limits<T>;
+    constexpr auto inf = limits::infinity();
+
+    ASSERT_TRUE(gko::isfinite(T{0}));
+    ASSERT_TRUE(gko::isfinite(-T{0}));
+    ASSERT_TRUE(gko::isfinite(T{1}));
+    ASSERT_FALSE(gko::isfinite(inf));
+    ASSERT_FALSE(gko::isfinite(-inf));
+    ASSERT_FALSE(gko::isfinite(limits::quiet_NaN()));
+    ASSERT_FALSE(gko::isfinite(limits::signaling_NaN()));
+    ASSERT_FALSE(gko::isfinite(inf - inf));    // results in nan
+    ASSERT_FALSE(gko::isfinite(inf / inf));    // results in nan
+    ASSERT_FALSE(gko::isfinite(inf * T{2}));   // results in inf
+    ASSERT_FALSE(gko::isfinite(T{1} / T{0}));  // results in inf
+    ASSERT_FALSE(gko::isfinite(T{0} / T{0}));  // results in nan
+}
+
+
+template <typename ComplexType>
+void test_complex_isfinite()
+{
+    static_assert(gko::is_complex_s<ComplexType>::value,
+                  "Template type must be a complex type.");
+    using T = gko::remove_complex<ComplexType>;
+    using c_type = ComplexType;
+    using limits = std::numeric_limits<T>;
+    constexpr auto inf = limits::infinity();
+    constexpr auto quiet_nan = limits::quiet_NaN();
+    constexpr auto signaling_nan = limits::signaling_NaN();
+
+    ASSERT_TRUE(gko::isfinite(c_type{T{0}, T{0}}));
+    ASSERT_TRUE(gko::isfinite(c_type{-T{0}, -T{0}}));
+    ASSERT_TRUE(gko::isfinite(c_type{T{1}, T{0}}));
+    ASSERT_TRUE(gko::isfinite(c_type{T{0}, T{1}}));
+    ASSERT_FALSE(gko::isfinite(c_type{inf, T{0}}));
+    ASSERT_FALSE(gko::isfinite(c_type{-inf, T{0}}));
+    ASSERT_FALSE(gko::isfinite(c_type{quiet_nan, T{0}}));
+    ASSERT_FALSE(gko::isfinite(c_type{signaling_nan, T{0}}));
+    ASSERT_FALSE(gko::isfinite(c_type{T{0}, inf}));
+    ASSERT_FALSE(gko::isfinite(c_type{T{0}, -inf}));
+    ASSERT_FALSE(gko::isfinite(c_type{T{0}, quiet_nan}));
+    ASSERT_FALSE(gko::isfinite(c_type{T{0}, signaling_nan}));
+}
+
+
+TEST(IsFinite, Float) { test_real_isfinite<float>(); }
+
+
+TEST(IsFinite, Double) { test_real_isfinite<double>(); }
+
+
+TEST(IsFinite, FloatComplex) { test_complex_isfinite<std::complex<float>>(); }
+
+
+TEST(IsFinite, DoubleComplex) { test_complex_isfinite<std::complex<double>>(); }
+
+
+}  // namespace

cuda/base/math.hpp

@@ -93,17 +93,125 @@ __device__ GKO_INLINE std::complex<double> one<std::complex<double>>()
 }
 
 
-// Since the `isfinite` function from CUDA is not in any namespace, it
-// will be used without the `using` keyword.
+// This first part is specific for clang and intel in combination with the nvcc
+// compiler from the toolkit older than 9.2.
+// Both want to use their `__builtin_isfinite` function, which is not present
+// as a __device__ function, so it results in a compiler error.
+// Here, `isfinite` is written by hand, which might not be as performant as the
+// intrinsic function from CUDA, but it compiles and works.
+#if defined(__CUDA_ARCH__) && defined(__CUDACC_VER_MAJOR__) &&     \
+    defined(__CUDACC_VER_MINOR__) &&                               \
+    (__CUDACC_VER_MAJOR__ * 1000 + __CUDACC_VER_MINOR__) < 9002 && \
+    (defined(__clang__) || defined(__ICC) || defined(__ICL))
 
-// For some CUDA versions, the complex variation for `isfinite` has to be
-// manually defined
-#if defined(GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE)
 
-GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE(thrust::complex<double>);
-GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE(thrust::complex<float>);
+namespace detail {
+
+
+/**
+ * This structure can be used to get the exponent mask of a given floating
+ * point type. Uses specialization to implement different types.
+ */
+template <typename T>
+struct mask_creator {};
+
+template <>
+struct mask_creator<float> {
+    using int_type = int32;
+    static constexpr int_type number_exponent_bits = 8;
+    static constexpr int_type number_significand_bits = 23;
+    // integer representation of a floating point number, where all exponent
+    // bits are set
+    static constexpr int_type exponent_mask =
+        ((int_type{1} << number_exponent_bits) - 1) << number_significand_bits;
+    static __device__ int_type reinterpret_int(const float &value)
+    {
+        return __float_as_int(value);
+    }
+};
+
+template <>
+struct mask_creator<double> {
+    using int_type = int64;
+    static constexpr int_type number_exponent_bits = 11;
+    static constexpr int_type number_significand_bits = 52;
+    // integer representation of a floating point number, where all exponent
+    // bits are set
+    static constexpr int_type exponent_mask =
+        ((int_type{1} << number_exponent_bits) - 1) << number_significand_bits;
+    static __device__ int_type reinterpret_int(const double &value)
+    {
+        return __double_as_longlong(value);
+    }
+};
+
+
+}  // namespace detail
+
 
-#endif  // defined(GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE)
+/**
+ * Checks if a given value is finite, meaning it is neither +/- infinity
+ * nor NaN.
+ *
+ * @internal  It checks if all exponent bits are set. If all are set, the
+ *            number either represents NaN or +/- infinity, meaning it is a
+ *            non-finite number.
+ *
+ * @param value  value to check
+ *
+ * returns `true` if the given value is finite, meaning it is neither
+ *         +/- infinity nor NaN.
+ */
+#define GKO_DEFINE_ISFINITE_FOR_TYPE(_type)                               \
+    GKO_INLINE __device__ bool isfinite(const _type &value)               \
+    {                                                                     \
+        constexpr auto mask = detail::mask_creator<_type>::exponent_mask; \
+        const auto re_int =                                               \
+            detail::mask_creator<_type>::reinterpret_int(value);          \
+        return (re_int & mask) != mask;                                   \
+    }
+
+GKO_DEFINE_ISFINITE_FOR_TYPE(float)
+GKO_DEFINE_ISFINITE_FOR_TYPE(double)
+#undef GKO_DEFINE_ISFINITE_FOR_TYPE
+
+
+/**
+ * Checks if all components of a complex value are finite, meaning they are
+ * neither +/- infinity nor NaN.
+ *
+ * @internal required for the clang compiler. This function will be used rather
+ *           than the `isfinite` function in the public `math.hpp` because
+ *           there is no template parameter, so it is prefered during lookup.
+ *
+ * @tparam T  complex type of the value to check
+ *
+ * @param value  complex value to check
+ *
+ * returns `true` if both components of the given value are finite, meaning
+ *         they are neither +/- infinity nor NaN.
+ */
+#define GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE(_type)              \
+    GKO_INLINE __device__ bool isfinite(const _type &value)      \
+    {                                                            \
+        return isfinite(value.real()) && isfinite(value.imag()); \
+    }
+
+GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE(thrust::complex<float>)
+GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE(thrust::complex<double>)
+#undef GKO_DEFINE_ISFINITE_FOR_COMPLEX_TYPE
+
+
+// For all other compiler in combination with CUDA, just use the provided
+// `isfinite` function
+#elif defined(__CUDA_ARCH__)
+
+
+// If it is compiled with the CUDA compiler, use their `isfinite`
+using ::isfinite;
+
+
+#endif  // defined(__CUDA_ARCH__)
 
 
 }  // namespace gko

cuda/test/base/CMakeLists.txt

@@ -1,2 +1,3 @@
 ginkgo_create_cuda_test(cuda_executor)
 ginkgo_create_cuda_test(exception_helpers)
+ginkgo_create_cuda_test(math)

cuda/test/base/math.cu

@@ -0,0 +1,145 @@
+/*******************************<GINKGO LICENSE>******************************
+Copyright (c) 2017-2019, the Ginkgo authors
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+******************************<GINKGO LICENSE>*******************************/
+
+#include <ginkgo/core/base/math.hpp>
+
+
+#include <cmath>
+#include <complex>
+#include <memory>
+
+
+#include <gtest/gtest.h>
+
+
+#include <ginkgo/core/base/array.hpp>
+#include <ginkgo/core/base/executor.hpp>
+
+
+#include "cuda/base/math.hpp"
+#include "cuda/base/types.hpp"
+
+
+namespace {
+
+
+template <typename T>
+__global__ void test_real_isfinite(bool *result)
+{
+    constexpr T inf = INFINITY;
+    bool test_true{};
+    bool test_false{};
+
+    test_true =
+        gko::isfinite(T{0}) && gko::isfinite(-T{0}) && gko::isfinite(T{1});
+    test_false = gko::isfinite(inf) || gko::isfinite(-inf) ||
+                 gko::isfinite(NAN) || gko::isfinite(inf - inf) ||
+                 gko::isfinite(inf / inf) || gko::isfinite(inf * T{2}) ||
+                 gko::isfinite(T{1} / T{0}) || gko::isfinite(T{0} / T{0});
+    *result = test_true && !test_false;
+}
+
+
+template <typename ComplexType>
+__global__ void test_complex_isfinite(bool *result)
+{
+    static_assert(gko::is_complex_s<ComplexType>::value,
+                  "Template type must be a complex type.");
+    using T = gko::remove_complex<ComplexType>;
+    using c_type = gko::kernels::cuda::cuda_type<ComplexType>;
+    constexpr T inf = INFINITY;
+    constexpr T quiet_nan = NAN;
+    bool test_true{};
+    bool test_false{};
+
+    test_true = gko::isfinite(c_type{T{0}, T{0}}) &&
+                gko::isfinite(c_type{-T{0}, -T{0}}) &&
+                gko::isfinite(c_type{T{1}, T{0}}) &&
+                gko::isfinite(c_type{T{0}, T{1}});
+    test_false =
+        gko::isfinite(c_type{inf, T{0}}) || gko::isfinite(c_type{-inf, T{0}}) ||
+        gko::isfinite(c_type{quiet_nan, T{0}}) ||
+        gko::isfinite(c_type{T{0}, inf}) || gko::isfinite(c_type{T{0}, -inf}) ||
+        gko::isfinite(c_type{T{0}, quiet_nan});
+    *result = test_true && !test_false;
+}
+
+
+class IsFinite : public ::testing::Test {
+protected:
+    IsFinite()
+        : ref(gko::ReferenceExecutor::create()),
+          cuda(gko::CudaExecutor::create(0, ref))
+    {}
+
+    template <typename T>
+    bool test_real_isfinite_kernel()
+    {
+        gko::Array<bool> result(cuda, 1);
+        test_real_isfinite<T><<<1, 1>>>(result.get_data());
+        result.set_executor(ref);
+        return *result.get_data();
+    }
+
+    template <typename T>
+    bool test_complex_isfinite_kernel()
+    {
+        gko::Array<bool> result(cuda, 1);
+        test_complex_isfinite<T><<<1, 1>>>(result.get_data());
+        result.set_executor(ref);
+        return *result.get_data();
+    }
+
+    std::shared_ptr<gko::ReferenceExecutor> ref;
+    std::shared_ptr<gko::CudaExecutor> cuda;
+};
+
+
+TEST_F(IsFinite, Float) { ASSERT_TRUE(test_real_isfinite_kernel<float>()); }
+
+
+TEST_F(IsFinite, Double) { ASSERT_TRUE(test_real_isfinite_kernel<double>()); }
+
+
+TEST_F(IsFinite, FloatComplex)
+{
+    ASSERT_TRUE(test_complex_isfinite_kernel<thrust::complex<float>>());
+}
+
+
+TEST_F(IsFinite, DoubleComplex)
+{
+    ASSERT_TRUE(test_complex_isfinite_kernel<thrust::complex<double>>());
+}
+
+
+}  // namespace