The first draft of simulated double.

This can be compiled with cuda 12.6.
It mainly focuses on the interface and forward the operation internally to double precision.
If the corresponding instruction is also removed by the vendor, some of them might be able to simulate by more operations (or allowing reinterpret), but some of them like atomic on 64 bits or memory control will reuqire hardware.

Author: yhmtsai

accessor/cuda_helper.hpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -25,6 +25,8 @@ namespace gko {
 
 class half;
 
+class custom_double;
+
 
 namespace acc {
 namespace detail {
@@ -40,6 +42,11 @@ struct cuda_type<gko::half> {
     using type = __half;
 };
 
+template <>
+struct cuda_type<double> {
+    using type = gko::custom_double;
+};
+
 // Unpack cv and reference / pointer qualifiers
 template <typename T>
 struct cuda_type<const T> {

accessor/reduced_row_major_reference.hpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -74,7 +74,7 @@ class reduced_storage
     operator=(arithmetic_type val) &&
     {
         storage_type* const GKO_ACC_RESTRICT r_ptr = ptr_;
-        *r_ptr = val;
+        *r_ptr = detail::implicit_explicit_conversion<storage_type>(val);
         return val;
     }
 

accessor/scaled_reduced_row_major_reference.hpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -76,7 +76,8 @@ class scaled_reduced_storage
     operator=(arithmetic_type val) &&
     {
         storage_type* const GKO_ACC_RESTRICT r_ptr = ptr_;
-        *r_ptr = val / scalar_;
+        *r_ptr =
+            detail::implicit_explicit_conversion<storage_type>(val / scalar_);
         return val;
     }
 

common/cuda_hip/base/device_matrix_data_kernels.cpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -27,7 +27,7 @@ namespace components {
 // Although gko::is_nonzero is constexpr, it still shows calling __device__ in
 // __host__
 template <typename T>
-GKO_INLINE __device__ constexpr bool is_nonzero(T value)
+GKO_INLINE __device__ constexpr bool is_nonzero_(T value)
 {
     return value != zero<T>();
 }
@@ -43,7 +43,7 @@ void remove_zeros(std::shared_ptr<const DefaultExecutor> exec,
     // count nonzeros
     auto nnz = thrust::count_if(
         thrust_policy(exec), value_ptr, value_ptr + size,
-        [] __device__(device_value_type value) { return is_nonzero(value); });
+        [] __device__(device_value_type value) { return is_nonzero_(value); });
     if (nnz < size) {
         using tuple_type =
             thrust::tuple<IndexType, IndexType, device_value_type>;
@@ -59,7 +59,7 @@ void remove_zeros(std::shared_ptr<const DefaultExecutor> exec,
                                as_device_type(new_values.get_data())));
         thrust::copy_if(thrust_policy(exec), it, it + size, out_it,
                         [] __device__(tuple_type entry) {
-                            return is_nonzero(thrust::get<2>(entry));
+                            return is_nonzero_(thrust::get<2>(entry));
                         });
         // swap out storage
         values = std::move(new_values);

common/cuda_hip/base/math.hpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -27,7 +27,7 @@
 
 
 #include "common/cuda_hip/base/thrust_macro.hpp"
-
+#include "core/base/custom_double.hpp"
 
 namespace gko {
 
@@ -67,6 +67,14 @@ struct device_numeric_limits<__half> {
     }
 };
 
+template <>
+GKO_INLINE constexpr custom_double one<custom_double>()
+{
+    constexpr auto bits = static_cast<uint64>(
+        0b0'01111111111'0000000000000000000000000000000000000000000000000000ull);
+    return custom_double::create_from_bits(bits);
+}
+
 
 namespace detail {
 
@@ -95,6 +103,9 @@ struct is_complex_impl<thrust::complex<T>> : public std::true_type {};
 template <>
 struct is_complex_or_scalar_impl<__half> : public std::true_type {};
 
+template <>
+struct is_complex_or_scalar_impl<gko::custom_double> : public std::true_type {};
+
 template <typename T>
 struct is_complex_or_scalar_impl<thrust::complex<T>>
     : public is_complex_or_scalar_impl<T> {};
@@ -125,6 +136,30 @@ GKO_ATTRIBUTES GKO_INLINE __half abs<__half>(const complex<__half>& z)
 }
 
 
+template <>
+GKO_ATTRIBUTES GKO_INLINE complex<gko::custom_double> sqrt<gko::custom_double>(
+    const complex<gko::custom_double>& a)
+{
+    auto result =
+        sqrt(complex<double>(gko::custom_double::custom_to_native(a.real()),
+                             gko::custom_double::custom_to_native(a.imag())));
+    return complex<gko::custom_double>(
+        gko::custom_double::to_custom(result.real()),
+        gko::custom_double::to_custom(result.imag()));
+}
+
+
+template <>
+GKO_ATTRIBUTES GKO_INLINE gko::custom_double abs<gko::custom_double>(
+    const complex<gko::custom_double>& z)
+{
+    auto result =
+        abs(complex<double>(gko::custom_double::custom_to_native(z.real()),
+                            gko::custom_double::custom_to_native(z.imag())));
+    return gko::custom_double::to_custom(result);
+}
+
+
 }  // namespace thrust
 GKO_THRUST_NAEMSPACE_POSTFIX
 
@@ -183,6 +218,33 @@ __device__ __forceinline__ bool is_finite(const thrust::complex<__half>& value)
     return is_finite(value.real()) && is_finite(value.imag());
 }
 
+
+__device__ __forceinline__ bool is_nan(const gko::custom_double& val)
+{
+    return is_nan(gko::custom_double::custom_to_native(val));
+}
+
+__device__ __forceinline__ bool is_nan(
+    const thrust::complex<gko::custom_double>& val)
+{
+    return is_nan(val.real()) || is_nan(val.imag());
+}
+
+__device__ __forceinline__ gko::custom_double abs(const gko::custom_double& val)
+{
+    return custom_double::to_custom(abs(custom_double::custom_to_native(val)));
+}
+
+__device__ __forceinline__ gko::custom_double sqrt(
+    const gko::custom_double& val)
+{
+    return custom_double::to_custom(sqrt(custom_double::custom_to_native(val)));
+}
+
+__device__ __forceinline__ bool is_finite(const gko::custom_double& value)
+{
+    return is_finite(custom_double::custom_to_native(value));
+}
 #endif
 
 

common/cuda_hip/components/atomic.hpp

@@ -1,4 +1,4 @@
-// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors
 //
 // SPDX-License-Identifier: BSD-3-Clause
 
@@ -240,6 +240,18 @@ __forceinline__ __device__ thrust::complex<double> atomic_add(
 }
 
 
+__forceinline__ __device__ thrust::complex<gko::custom_double> atomic_add(
+    thrust::complex<gko::custom_double>* __restrict__ address,
+    thrust::complex<gko::custom_double> val)
+{
+    auto addr = reinterpret_cast<gko::custom_double*>(address);
+    // Separate to real part and imag part
+    auto real = atomic_add(addr, val.real());
+    auto imag = atomic_add(addr + 1, val.imag());
+    return {real, imag};
+}
+
+
 }  // namespace GKO_DEVICE_NAMESPACE
 }  // namespace kernels
 }  // namespace gko