The draft of simulated double.

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