1888/precision handling

core/CMakeLists.txt

@@ -37,10 +37,12 @@ target_sources(
         base/device_matrix_data.cpp
         base/executor.cpp
         base/index_set.cpp
+        base/lin_op.cpp
         base/memory.cpp
         base/mpi.cpp
         base/mtx_io.cpp
         base/perturbation.cpp
+        base/precision.cpp
         base/segmented_array.cpp
         base/timer.cpp
         base/version.cpp

core/base/block_operator.cpp

@@ -153,14 +153,14 @@ std::unique_ptr<BlockOperator> BlockOperator::create(
 
 
 BlockOperator::BlockOperator(std::shared_ptr<const Executor> exec)
-    : LinOp(std::move(exec))
+    : LinOp(std::move(exec), dim<2>{}, precision::any)
 {}
 
 
 BlockOperator::BlockOperator(
     std::shared_ptr<const Executor> exec,
     std::vector<std::vector<std::shared_ptr<const LinOp>>> blocks)
-    : LinOp(exec, compute_global_size(blocks)),
+    : LinOp(exec, compute_global_size(blocks), precision::any),
       block_size_(blocks.empty()
                       ? dim<2>{}
                       : dim<2>(blocks.size(), blocks.front().size())),
@@ -246,14 +246,14 @@ void BlockOperator::apply_impl(const LinOp* alpha, const LinOp* b,
 
 
 BlockOperator::BlockOperator(const BlockOperator& other)
-    : LinOp(other.get_executor())
+    : LinOp(other.get_executor(), dim<2>{}, precision::any)
 {
     *this = other;
 }
 
 
 BlockOperator::BlockOperator(BlockOperator&& other) noexcept
-    : LinOp(other.get_executor())
+    : LinOp(other.get_executor(), dim<2>{}, precision::any)
 {
     *this = std::move(other);
 }

core/base/lin_op.cpp

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: 2026 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <ginkgo/core/base/lin_op.hpp>
+
+namespace gko {
+
+
+LinOp::LinOp(LinOp&& other)
+    : PolymorphicObject(std::move(other)),
+      size_{std::exchange(other.size_, dim<2>{})},
+      value_t_(other.value_t_)
+{}
+
+
+LinOp::LinOp(std::shared_ptr<const Executor> exec, const dim<2>& size,
+             precision p)
+    : PolymorphicObject(exec), size_{size}, value_t_(p)
+{}
+
+
+precision LinOp::get_precision() const noexcept { return value_t_; }
+
+
+void LinOp::set_precision(precision p) noexcept { value_t_ = p; }
+
+
+}  // namespace gko

core/base/perturbation.cpp

@@ -68,7 +68,7 @@ Perturbation<ValueType>::Perturbation(Perturbation&& other)
 
 template <typename ValueType>
 Perturbation<ValueType>::Perturbation(std::shared_ptr<const Executor> exec)
-    : LinOp(std::move(exec))
+    : LinOp(std::move(exec), dim<2>{}, precision_v<ValueType>)
 {}
 
 
@@ -87,7 +87,8 @@ template <typename ValueType>
 Perturbation<ValueType>::Perturbation(std::shared_ptr<const LinOp> scalar,
                                       std::shared_ptr<const LinOp> basis,
                                       std::shared_ptr<const LinOp> projector)
-    : LinOp(basis->get_executor(), gko::dim<2>{basis->get_size()[0]}),
+    : LinOp(basis->get_executor(), gko::dim<2>{basis->get_size()[0]},
+            precision_v<ValueType>),
       basis_{std::move(basis)},
       projector_{std::move(projector)},
       scalar_{std::move(scalar)}

core/base/precision.cpp

@@ -0,0 +1,131 @@
+// SPDX-FileCopyrightText: 2026 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include "ginkgo/core/base/precision.hpp"
+
+#include <ginkgo/core/base/exception_helpers.hpp>
+
+
+namespace gko {
+
+
+precision as_real(precision p)
+{
+    if (is_real(p)) {
+        return p;
+    }
+    switch (p) {
+    case precision::complex_fp32:
+        return precision::fp32;
+    case precision::complex_fp64:
+        return precision::fp64;
+#if GINKGO_ENABLE_HALF
+    case precision::complex_fp16:
+        return precision::fp16;
+#endif
+#if GINKGO_ENABLE_BFLOAT16
+    case precision::complex_bf16:
+        return precision::bf16;
+#endif
+    default:
+        GKO_INVALID_STATE("Unsupported precision");
+    }
+}
+
+
+precision as_complex(precision p)
+{
+    if (is_complex(p)) {
+        return p;
+    }
+    switch (p) {
+    case precision::fp32:
+        return precision::complex_fp32;
+    case precision::fp64:
+        return precision::complex_fp64;
+#if GINKGO_ENABLE_HALF
+    case precision::fp16:
+        return precision::complex_fp16;
+#endif
+#if GINKGO_ENABLE_BFLOAT16
+    case precision::bf16:
+        return precision::complex_bf16;
+#endif
+    default:
+        GKO_INVALID_STATE("Unsupported precision");
+    }
+}
+
+
+std::string to_string(precision p)
+{
+    switch (p) {
+    case precision::fp32:
+        return "fp32";
+    case precision::complex_fp32:
+        return "complex_fp32";
+    case precision::fp64:
+        return "fp64";
+    case precision::complex_fp64:
+        return "complex_fp64";
+#if GINKGO_ENABLE_HALF
+    case precision::fp16:
+        return "fp16";
+    case precision::complex_fp16:
+        return "complex_fp16";
+#endif
+#if GINKGO_ENABLE_BFLOAT16
+    case precision::bf16:
+        return "bf16";
+    case precision::complex_bf16:
+        return "complex_bf16";
+#endif
+    case precision::any:
+        return "any";
+    case precision::none:
+        return "none";
+    default:
+        GKO_INVALID_STATE("Unsupported precision");
+    }
+}
+
+
+std::variant<
+#if GINKGO_ENABLE_HALF
+    half, std::complex<half>,
+#endif
+#if GINKGO_ENABLE_BFLOAT16
+    bfloat16, std::complex<bfloat16>,
+#endif
+    float, std::complex<float>, double, std::complex<double>>
+precision_to_variant(precision p)
+{
+    switch (p) {
+#if GINKGO_ENABLE_HALF
+    case precision::fp16:
+        return half{};
+    case precision::complex_fp16:
+        return std::complex<half>{};
+#endif
+#if GINKGO_ENABLE_BFLOAT16
+    case precision::bf16:
+        return bfloat16{};
+    case precision::complex_bf16:
+        return std::complex<bfloat16>{};
+#endif
+    case precision::fp32:
+        return float{};
+    case precision::complex_fp32:
+        return std::complex<float>{};
+    case precision::fp64:
+        return double{};
+    case precision::complex_fp64:
+        return std::complex<double>{};
+    default:
+        GKO_INVALID_STATE("Unsupported precision");
+    }
+}
+
+
+}  // namespace gko

core/distributed/matrix.cpp

@@ -89,7 +89,7 @@ Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(
     std::shared_ptr<const RowGatherer<LocalIndexType>> row_gather_template,
     ptr_param<const LinOp> diag_matrix_template,
     ptr_param<const LinOp> off_diag_matrix_template)
-    : LinOp{exec},
+    : LinOp{exec, dim<2>{}, precision_v<ValueType>},
       DistributedBase{row_gather_template->get_communicator()},
       row_gatherer_{clone(exec, row_gather_template)},
       imap_{exec},
@@ -109,7 +109,7 @@ template <typename ValueType, typename LocalIndexType, typename GlobalIndexType>
 Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(
     std::shared_ptr<const Executor> exec, mpi::communicator comm, dim<2> size,
     std::shared_ptr<LinOp> diag_linop)
-    : LinOp{exec},
+    : LinOp{exec, dim<2>{}, precision_v<ValueType>},
       DistributedBase{comm},
       row_gatherer_{RowGatherer<LocalIndexType>::create(
           exec, mpi::detail::create_default_collective_communicator(comm))},
@@ -132,7 +132,7 @@ Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(
     std::shared_ptr<const Executor> exec, mpi::communicator comm,
     index_map<LocalIndexType, GlobalIndexType> imap,
     std::shared_ptr<LinOp> diag_linop, std::shared_ptr<LinOp> off_diag_linop)
-    : LinOp{exec},
+    : LinOp{exec, dim<2>{}, precision_v<ValueType>},
       DistributedBase{comm},
       row_gatherer_(RowGatherer<LocalIndexType>::create(
           exec,
@@ -748,7 +748,7 @@ void Matrix<ValueType, LocalIndexType, GlobalIndexType>::row_scale(
 
 template <typename ValueType, typename LocalIndexType, typename GlobalIndexType>
 Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(const Matrix& other)
-    : LinOp{other.get_executor()},
+    : LinOp{other.get_executor(), dim<2>{}, precision_v<ValueType>},
       DistributedBase{other.get_communicator()},
       row_gatherer_{RowGatherer<LocalIndexType>::create(
           other.get_executor(), other.get_communicator())},
@@ -762,7 +762,7 @@ Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(const Matrix& other)
 template <typename ValueType, typename LocalIndexType, typename GlobalIndexType>
 Matrix<ValueType, LocalIndexType, GlobalIndexType>::Matrix(
     Matrix&& other) noexcept
-    : LinOp{other.get_executor()},
+    : LinOp{other.get_executor(), dim<2>{}, precision_v<ValueType>},
       DistributedBase{other.get_communicator()},
       row_gatherer_{RowGatherer<LocalIndexType>::create(
           other.get_executor(), other.get_communicator())},

core/factorization/factorization.cpp

@@ -243,7 +243,7 @@ Factorization<ValueType, IndexType>::operator=(Factorization&& fact)
 template <typename ValueType, typename IndexType>
 Factorization<ValueType, IndexType>::Factorization(
     std::shared_ptr<const Executor> exec)
-    : LinOp{exec},
+    : LinOp{exec, dim<2>{}, precision_v<ValueType>},
       storage_type_{storage_type::empty},
       factors_{Composition<ValueType>::create(exec)}
 {}
@@ -252,7 +252,8 @@ Factorization<ValueType, IndexType>::Factorization(
 template <typename ValueType, typename IndexType>
 Factorization<ValueType, IndexType>::Factorization(
     std::unique_ptr<Composition<ValueType>> factors, storage_type type)
-    : LinOp{factors->get_executor(), factors->get_size()},
+    : LinOp{factors->get_executor(), factors->get_size(),
+            precision_v<ValueType>},
       storage_type_{type},
       factors_{std::move(factors)}
 {}

core/matrix/coo.cpp

@@ -91,7 +91,7 @@ Coo<ValueType, IndexType>::create_const(
 template <typename ValueType, typename IndexType>
 Coo<ValueType, IndexType>::Coo(std::shared_ptr<const Executor> exec,
                                const dim<2>& size, size_type num_nonzeros)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       values_(exec, num_nonzeros),
       col_idxs_(exec, num_nonzeros),
       row_idxs_(exec, num_nonzeros)
@@ -103,7 +103,7 @@ Coo<ValueType, IndexType>::Coo(std::shared_ptr<const Executor> exec,
                                const dim<2>& size, array<value_type> values,
                                array<index_type> col_idxs,
                                array<index_type> row_idxs)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       values_{exec, std::move(values)},
       col_idxs_{exec, std::move(col_idxs)},
       row_idxs_{exec, std::move(row_idxs)}

core/matrix/csr.cpp

@@ -278,7 +278,7 @@ template <typename ValueType, typename IndexType>
 Csr<ValueType, IndexType>::Csr(std::shared_ptr<const Executor> exec,
                                const dim<2>& size, size_type num_nonzeros,
                                csr::spmv_strategy strategy)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       strategy_(strategy),
       values_(exec, num_nonzeros),
       col_idxs_(exec, num_nonzeros),
@@ -296,7 +296,7 @@ Csr<ValueType, IndexType>::Csr(std::shared_ptr<const Executor> exec,
                                array<index_type> col_idxs,
                                array<index_type> row_ptrs,
                                csr::spmv_strategy strategy)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       strategy_(strategy),
       values_{exec, std::move(values)},
       col_idxs_{exec, std::move(col_idxs)},

core/matrix/dense.cpp

@@ -936,7 +936,7 @@ Dense<ValueType>& Dense<ValueType>::operator=(Dense&& other)
 template <typename ValueType>
 Dense<ValueType>::Dense(std::shared_ptr<const Executor> exec,
                         const dim<2>& size, size_type stride)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       stride_(stride == 0 ? size[1] : stride),
       values_(exec, size[0] * stride_)
 {}
@@ -946,7 +946,7 @@ template <typename ValueType>
 Dense<ValueType>::Dense(std::shared_ptr<const Executor> exec,
                         const dim<2>& size, array<value_type> values,
                         size_type stride)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       stride_(stride == 0 ? size[1] : stride),
       values_(exec, std::move(values))
 {

core/matrix/diagonal.cpp

@@ -368,14 +368,15 @@ Diagonal<ValueType>::compute_absolute() const
 template <typename ValueType>
 Diagonal<ValueType>::Diagonal(std::shared_ptr<const Executor> exec,
                               size_type size)
-    : LinOp(exec, dim<2>{size}), values_(exec, size)
+    : LinOp(exec, dim<2>{size}, precision_v<ValueType>), values_(exec, size)
 {}
 
 
 template <typename ValueType>
 Diagonal<ValueType>::Diagonal(std::shared_ptr<const Executor> exec,
                               const size_type size, array<value_type> values)
-    : LinOp(exec, dim<2>(size)), values_{exec, std::move(values)}
+    : LinOp(exec, dim<2>(size), precision_v<ValueType>),
+      values_{exec, std::move(values)}
 {
     GKO_ENSURE_COMPATIBLE_BOUNDS(size, values_.get_size());
 }

core/matrix/ell.cpp

@@ -449,7 +449,7 @@ Ell<ValueType, IndexType>::Ell(std::shared_ptr<const Executor> exec,
                                const dim<2>& size,
                                size_type num_stored_elements_per_row,
                                size_type stride)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       num_stored_elements_per_row_(num_stored_elements_per_row),
       stride_(stride == 0 ? size[0] : stride),
       values_(exec, stride_ * num_stored_elements_per_row),
@@ -463,7 +463,7 @@ Ell<ValueType, IndexType>::Ell(std::shared_ptr<const Executor> exec,
                                array<index_type> col_idxs,
                                size_type num_stored_elements_per_row,
                                size_type stride)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       num_stored_elements_per_row_{num_stored_elements_per_row},
       stride_{stride},
       values_{exec, std::move(values)},

core/matrix/fbcsr.cpp

@@ -494,7 +494,7 @@ template <typename ValueType, typename IndexType>
 Fbcsr<ValueType, IndexType>::Fbcsr(std::shared_ptr<const Executor> exec,
                                    const dim<2>& size, size_type num_nonzeros,
                                    int block_size)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       bs_{block_size},
       values_(exec, num_nonzeros),
       col_idxs_(exec,
@@ -512,7 +512,7 @@ Fbcsr<ValueType, IndexType>::Fbcsr(std::shared_ptr<const Executor> exec,
                                    array<value_type> values,
                                    array<index_type> col_idxs,
                                    array<index_type> row_ptrs)
-    : LinOp(exec, size),
+    : LinOp(exec, size, precision_v<ValueType>),
       bs_{block_size},
       values_{exec, std::move(values)},
       col_idxs_{exec, std::move(col_idxs)},