diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index a2eed9a48d6..40ef99637da 100644 --- a/core/CMakeLists.txt +++ b/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 diff --git a/core/base/block_operator.cpp b/core/base/block_operator.cpp index 700513ad9c2..8bb0a55a0d9 100644 --- a/core/base/block_operator.cpp +++ b/core/base/block_operator.cpp @@ -153,14 +153,14 @@ std::unique_ptr BlockOperator::create( BlockOperator::BlockOperator(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision::any) {} BlockOperator::BlockOperator( std::shared_ptr exec, std::vector>> 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); } diff --git a/core/base/lin_op.cpp b/core/base/lin_op.cpp new file mode 100644 index 00000000000..3db49ce32d1 --- /dev/null +++ b/core/base/lin_op.cpp @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +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 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 diff --git a/core/base/perturbation.cpp b/core/base/perturbation.cpp index 5b48ba0ddb1..45bfa8f521a 100644 --- a/core/base/perturbation.cpp +++ b/core/base/perturbation.cpp @@ -68,7 +68,7 @@ Perturbation::Perturbation(Perturbation&& other) template Perturbation::Perturbation(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} @@ -87,7 +87,8 @@ template Perturbation::Perturbation(std::shared_ptr scalar, std::shared_ptr basis, std::shared_ptr 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), basis_{std::move(basis)}, projector_{std::move(projector)}, scalar_{std::move(scalar)} diff --git a/core/base/precision.cpp b/core/base/precision.cpp new file mode 100644 index 00000000000..a81594cb5a2 --- /dev/null +++ b/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 + + +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, +#endif +#if GINKGO_ENABLE_BFLOAT16 + bfloat16, std::complex, +#endif + float, std::complex, double, std::complex> +precision_to_variant(precision p) +{ + switch (p) { +#if GINKGO_ENABLE_HALF + case precision::fp16: + return half{}; + case precision::complex_fp16: + return std::complex{}; +#endif +#if GINKGO_ENABLE_BFLOAT16 + case precision::bf16: + return bfloat16{}; + case precision::complex_bf16: + return std::complex{}; +#endif + case precision::fp32: + return float{}; + case precision::complex_fp32: + return std::complex{}; + case precision::fp64: + return double{}; + case precision::complex_fp64: + return std::complex{}; + default: + GKO_INVALID_STATE("Unsupported precision"); + } +} + + +} // namespace gko diff --git a/core/distributed/matrix.cpp b/core/distributed/matrix.cpp index 60c0f4ca64b..fa191cf13b0 100644 --- a/core/distributed/matrix.cpp +++ b/core/distributed/matrix.cpp @@ -89,7 +89,7 @@ Matrix::Matrix( std::shared_ptr> row_gather_template, ptr_param diag_matrix_template, ptr_param off_diag_matrix_template) - : LinOp{exec}, + : LinOp{exec, dim<2>{}, precision_v}, DistributedBase{row_gather_template->get_communicator()}, row_gatherer_{clone(exec, row_gather_template)}, imap_{exec}, @@ -109,7 +109,7 @@ template Matrix::Matrix( std::shared_ptr exec, mpi::communicator comm, dim<2> size, std::shared_ptr diag_linop) - : LinOp{exec}, + : LinOp{exec, dim<2>{}, precision_v}, DistributedBase{comm}, row_gatherer_{RowGatherer::create( exec, mpi::detail::create_default_collective_communicator(comm))}, @@ -132,7 +132,7 @@ Matrix::Matrix( std::shared_ptr exec, mpi::communicator comm, index_map imap, std::shared_ptr diag_linop, std::shared_ptr off_diag_linop) - : LinOp{exec}, + : LinOp{exec, dim<2>{}, precision_v}, DistributedBase{comm}, row_gatherer_(RowGatherer::create( exec, @@ -748,7 +748,7 @@ void Matrix::row_scale( template Matrix::Matrix(const Matrix& other) - : LinOp{other.get_executor()}, + : LinOp{other.get_executor(), dim<2>{}, precision_v}, DistributedBase{other.get_communicator()}, row_gatherer_{RowGatherer::create( other.get_executor(), other.get_communicator())}, @@ -762,7 +762,7 @@ Matrix::Matrix(const Matrix& other) template Matrix::Matrix( Matrix&& other) noexcept - : LinOp{other.get_executor()}, + : LinOp{other.get_executor(), dim<2>{}, precision_v}, DistributedBase{other.get_communicator()}, row_gatherer_{RowGatherer::create( other.get_executor(), other.get_communicator())}, diff --git a/core/factorization/factorization.cpp b/core/factorization/factorization.cpp index 3e1d8c7b077..d330f9cdfc2 100644 --- a/core/factorization/factorization.cpp +++ b/core/factorization/factorization.cpp @@ -243,7 +243,7 @@ Factorization::operator=(Factorization&& fact) template Factorization::Factorization( std::shared_ptr exec) - : LinOp{exec}, + : LinOp{exec, dim<2>{}, precision_v}, storage_type_{storage_type::empty}, factors_{Composition::create(exec)} {} @@ -252,7 +252,8 @@ Factorization::Factorization( template Factorization::Factorization( std::unique_ptr> factors, storage_type type) - : LinOp{factors->get_executor(), factors->get_size()}, + : LinOp{factors->get_executor(), factors->get_size(), + precision_v}, storage_type_{type}, factors_{std::move(factors)} {} diff --git a/core/matrix/coo.cpp b/core/matrix/coo.cpp index 4eb0067d740..eb8d0656581 100644 --- a/core/matrix/coo.cpp +++ b/core/matrix/coo.cpp @@ -91,7 +91,7 @@ Coo::create_const( template Coo::Coo(std::shared_ptr exec, const dim<2>& size, size_type num_nonzeros) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), values_(exec, num_nonzeros), col_idxs_(exec, num_nonzeros), row_idxs_(exec, num_nonzeros) @@ -103,7 +103,7 @@ Coo::Coo(std::shared_ptr exec, const dim<2>& size, array values, array col_idxs, array row_idxs) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), values_{exec, std::move(values)}, col_idxs_{exec, std::move(col_idxs)}, row_idxs_{exec, std::move(row_idxs)} diff --git a/core/matrix/csr.cpp b/core/matrix/csr.cpp index e70a3739692..a40b8340de2 100644 --- a/core/matrix/csr.cpp +++ b/core/matrix/csr.cpp @@ -278,7 +278,7 @@ template Csr::Csr(std::shared_ptr exec, const dim<2>& size, size_type num_nonzeros, csr::spmv_strategy strategy) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), strategy_(strategy), values_(exec, num_nonzeros), col_idxs_(exec, num_nonzeros), @@ -296,7 +296,7 @@ Csr::Csr(std::shared_ptr exec, array col_idxs, array row_ptrs, csr::spmv_strategy strategy) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), strategy_(strategy), values_{exec, std::move(values)}, col_idxs_{exec, std::move(col_idxs)}, diff --git a/core/matrix/dense.cpp b/core/matrix/dense.cpp index 9b3a08a0b25..bcb3bdbeecf 100644 --- a/core/matrix/dense.cpp +++ b/core/matrix/dense.cpp @@ -936,7 +936,7 @@ Dense& Dense::operator=(Dense&& other) template Dense::Dense(std::shared_ptr exec, const dim<2>& size, size_type stride) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), stride_(stride == 0 ? size[1] : stride), values_(exec, size[0] * stride_) {} @@ -946,7 +946,7 @@ template Dense::Dense(std::shared_ptr exec, const dim<2>& size, array values, size_type stride) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), stride_(stride == 0 ? size[1] : stride), values_(exec, std::move(values)) { diff --git a/core/matrix/diagonal.cpp b/core/matrix/diagonal.cpp index a24d1c42d37..7e1cffd224e 100644 --- a/core/matrix/diagonal.cpp +++ b/core/matrix/diagonal.cpp @@ -368,14 +368,15 @@ Diagonal::compute_absolute() const template Diagonal::Diagonal(std::shared_ptr exec, size_type size) - : LinOp(exec, dim<2>{size}), values_(exec, size) + : LinOp(exec, dim<2>{size}, precision_v), values_(exec, size) {} template Diagonal::Diagonal(std::shared_ptr exec, const size_type size, array values) - : LinOp(exec, dim<2>(size)), values_{exec, std::move(values)} + : LinOp(exec, dim<2>(size), precision_v), + values_{exec, std::move(values)} { GKO_ENSURE_COMPATIBLE_BOUNDS(size, values_.get_size()); } diff --git a/core/matrix/ell.cpp b/core/matrix/ell.cpp index 5314bde47a2..1ca1a014c4e 100644 --- a/core/matrix/ell.cpp +++ b/core/matrix/ell.cpp @@ -449,7 +449,7 @@ Ell::Ell(std::shared_ptr exec, const dim<2>& size, size_type num_stored_elements_per_row, size_type stride) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), 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::Ell(std::shared_ptr exec, array col_idxs, size_type num_stored_elements_per_row, size_type stride) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), num_stored_elements_per_row_{num_stored_elements_per_row}, stride_{stride}, values_{exec, std::move(values)}, diff --git a/core/matrix/fbcsr.cpp b/core/matrix/fbcsr.cpp index b422ceec4d1..a28ef70224c 100644 --- a/core/matrix/fbcsr.cpp +++ b/core/matrix/fbcsr.cpp @@ -494,7 +494,7 @@ template Fbcsr::Fbcsr(std::shared_ptr exec, const dim<2>& size, size_type num_nonzeros, int block_size) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), bs_{block_size}, values_(exec, num_nonzeros), col_idxs_(exec, @@ -512,7 +512,7 @@ Fbcsr::Fbcsr(std::shared_ptr exec, array values, array col_idxs, array row_ptrs) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), bs_{block_size}, values_{exec, std::move(values)}, col_idxs_{exec, std::move(col_idxs)}, diff --git a/core/matrix/fft.cpp b/core/matrix/fft.cpp index 95dfee84428..beb5a1a12ab 100644 --- a/core/matrix/fft.cpp +++ b/core/matrix/fft.cpp @@ -187,7 +187,9 @@ void Fft::apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, Fft::Fft(std::shared_ptr exec, size_type size, bool inverse) - : LinOp(exec, dim<2>{size}), buffer_{exec}, inverse_{inverse} + : LinOp(exec, dim<2>{size}, precision::any), + buffer_{exec}, + inverse_{inverse} {} @@ -307,7 +309,7 @@ std::unique_ptr Fft2::create(std::shared_ptr exec, Fft2::Fft2(std::shared_ptr exec, size_type size1, size_type size2, bool inverse) - : LinOp(exec, dim<2>{size1 * size2}), + : LinOp(exec, dim<2>{size1 * size2}, precision::any), buffer_{exec}, fft_size_{size1, size2}, inverse_{inverse} @@ -421,7 +423,7 @@ std::unique_ptr Fft3::create(std::shared_ptr exec, Fft3::Fft3(std::shared_ptr exec, size_type size1, size_type size2, size_type size3, bool inverse) - : LinOp(exec, dim<2>{size1 * size2 * size3}), + : LinOp(exec, dim<2>{size1 * size2 * size3}, precision::any), buffer_{exec}, fft_size_{size1, size2, size3}, inverse_{inverse} diff --git a/core/matrix/hybrid.cpp b/core/matrix/hybrid.cpp index 98e03b755a2..b3314e425c1 100644 --- a/core/matrix/hybrid.cpp +++ b/core/matrix/hybrid.cpp @@ -107,7 +107,7 @@ Hybrid::Hybrid(std::shared_ptr exec, size_type num_stored_elements_per_row, size_type stride, size_type num_nonzeros, std::shared_ptr strategy) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), ell_(ell_type::create(exec, size, num_stored_elements_per_row, stride)), coo_(coo_type::create(exec, size, num_nonzeros)), strategy_(strategy ? std::move(strategy) : std::make_shared()) diff --git a/core/matrix/identity.cpp b/core/matrix/identity.cpp index aea7b119380..85ec47a9e3d 100644 --- a/core/matrix/identity.cpp +++ b/core/matrix/identity.cpp @@ -60,7 +60,7 @@ std::unique_ptr Identity::conj_transpose() const template Identity::Identity(std::shared_ptr exec, size_type size) - : LinOp(exec, dim<2>{size}) + : LinOp(exec, dim<2>{size}, precision_v) {} diff --git a/core/matrix/multivector.cpp b/core/matrix/multivector.cpp index 704f54d1522..a71c0a57227 100644 --- a/core/matrix/multivector.cpp +++ b/core/matrix/multivector.cpp @@ -1760,7 +1760,7 @@ std::unique_ptr> MultiVector::as_dense_view() template MultiVector::MultiVector(std::shared_ptr exec, const dim<2>& size, size_type stride) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), stride_(stride == 0 ? size[1] : stride), values_(exec, size[0] * stride_) {} @@ -1770,7 +1770,9 @@ template MultiVector::MultiVector(std::shared_ptr exec, const dim<2>& size, array values, size_type stride) - : LinOp(exec, size), stride_{stride}, values_{exec, std::move(values)} + : LinOp(exec, size, precision_v), + stride_{stride}, + values_{exec, std::move(values)} { if (size[0] > 0 && size[1] > 0) { GKO_ENSURE_IN_BOUNDS((size[0] - 1) * stride + size[1] - 1, diff --git a/core/matrix/permutation.cpp b/core/matrix/permutation.cpp index 5652a9cc943..d30f08cb7a4 100644 --- a/core/matrix/permutation.cpp +++ b/core/matrix/permutation.cpp @@ -184,14 +184,14 @@ Permutation::create_const( template Permutation::Permutation(std::shared_ptr exec, size_type size) - : LinOp(exec, dim<2>{size}), permutation_{exec, size} + : LinOp(exec, dim<2>{size}, precision::any), permutation_{exec, size} {} template Permutation::Permutation(std::shared_ptr exec, array permutation_indices) - : LinOp(exec, dim<2>{permutation_indices.get_size()}), + : LinOp(exec, dim<2>{permutation_indices.get_size()}, precision::any), permutation_{exec, std::move(permutation_indices)} {} diff --git a/core/matrix/row_gatherer.cpp b/core/matrix/row_gatherer.cpp index 14ae32f4d62..0ea350acde1 100644 --- a/core/matrix/row_gatherer.cpp +++ b/core/matrix/row_gatherer.cpp @@ -17,7 +17,7 @@ namespace matrix { template RowGatherer::RowGatherer(std::shared_ptr exec, const dim<2>& size) - : LinOp(exec, size), row_idxs_(exec, size[0]) + : LinOp(exec, size, precision::any), row_idxs_(exec, size[0]) {} @@ -25,7 +25,7 @@ template RowGatherer::RowGatherer(std::shared_ptr exec, const dim<2>& size, array row_idxs) - : LinOp(exec, size), row_idxs_{exec, std::move(row_idxs)} + : LinOp(exec, size, precision::any), row_idxs_{exec, std::move(row_idxs)} { GKO_ASSERT_EQ(size[0], row_idxs_.get_size()); } diff --git a/core/matrix/sellp.cpp b/core/matrix/sellp.cpp index 8377a0f09c2..7e0d5c3a7b2 100644 --- a/core/matrix/sellp.cpp +++ b/core/matrix/sellp.cpp @@ -117,7 +117,7 @@ Sellp::Sellp(std::shared_ptr exec, const dim<2>& size, size_type slice_size, size_type stride_factor, size_type total_cols) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), values_(exec, slice_size * total_cols), col_idxs_(exec, slice_size * total_cols), slice_lengths_(exec, ceildiv(size[0], slice_size)), diff --git a/core/matrix/sparsity_csr.cpp b/core/matrix/sparsity_csr.cpp index 403848c3f4e..f5551e7847c 100644 --- a/core/matrix/sparsity_csr.cpp +++ b/core/matrix/sparsity_csr.cpp @@ -133,7 +133,7 @@ template SparsityCsr::SparsityCsr( std::shared_ptr exec, const dim<2>& size, size_type num_nonzeros) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), col_idxs_(exec, num_nonzeros), row_ptrs_(exec, size[0] + 1), value_(exec, {one()}) @@ -146,7 +146,7 @@ template SparsityCsr::SparsityCsr( std::shared_ptr exec, const dim<2>& size, array col_idxs, array row_ptrs, value_type value) - : LinOp(exec, size), + : LinOp(exec, size, precision_v), col_idxs_{exec, std::move(col_idxs)}, row_ptrs_{exec, std::move(row_ptrs)}, value_{exec, {value}} @@ -158,7 +158,7 @@ SparsityCsr::SparsityCsr( template SparsityCsr::SparsityCsr( std::shared_ptr exec, std::shared_ptr matrix) - : LinOp(exec, matrix->get_size()) + : LinOp(exec, matrix->get_size(), precision_v) { auto tmp_ = copy_and_convert_to(exec, matrix); this->copy_from(tmp_); diff --git a/core/preconditioner/ic.cpp b/core/preconditioner/ic.cpp index d12cb783e59..4c8e3d1897d 100644 --- a/core/preconditioner/ic.cpp +++ b/core/preconditioner/ic.cpp @@ -160,14 +160,15 @@ void Ic::apply_impl(const LinOp* alpha, const LinOp* b, template Ic::Ic(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} template Ic::Ic(const Factory* factory, std::shared_ptr lin_op) - : LinOp(factory->get_executor(), lin_op->get_size()), + : LinOp(factory->get_executor(), lin_op->get_size(), + precision_v), parameters_{factory->get_parameters()} { auto comp = diff --git a/core/preconditioner/ilu.cpp b/core/preconditioner/ilu.cpp index 1d560d3a86e..8eb0fadfe25 100644 --- a/core/preconditioner/ilu.cpp +++ b/core/preconditioner/ilu.cpp @@ -207,14 +207,15 @@ void Ilu::apply_impl(const LinOp* alpha, template Ilu::Ilu( std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} template Ilu::Ilu( const Factory* factory, std::shared_ptr lin_op) - : LinOp(factory->get_executor(), lin_op->get_size()), + : LinOp(factory->get_executor(), lin_op->get_size(), + precision_v), parameters_{factory->get_parameters()} { auto comp = diff --git a/core/preconditioner/isai.cpp b/core/preconditioner/isai.cpp index 4d8e57bbd92..6a8fa9bd8fa 100644 --- a/core/preconditioner/isai.cpp +++ b/core/preconditioner/isai.cpp @@ -120,6 +120,49 @@ Isai::parse( return params; } + +template +Isai::Isai( + const Factory* factory, std::shared_ptr system_matrix) + + : LinOp(factory->get_executor(), system_matrix->get_size(), + precision_v), + parameters_{factory->get_parameters()} +{ + const auto skip_sorting = parameters_.skip_sorting; + const auto power = parameters_.sparsity_power; + const auto excess_limit = parameters_.excess_limit; + generate_inverse(system_matrix, skip_sorting, power, excess_limit, + static_cast>( + parameters_.excess_solver_reduction)); + if (IsaiType == isai_type::spd) { + auto inv = share(as(approximate_inverse_)); + auto inv_transp = share(inv->conj_transpose()); + approximate_inverse_ = Composition::create(inv_transp, inv); + } +} + + +template +void Isai::apply_impl(const LinOp* b, + LinOp* x) const + +{ + approximate_inverse_->apply(b, x); +} + + +template +void Isai::apply_impl(const LinOp* alpha, + const LinOp* b, + const LinOp* beta, + LinOp* x) const + +{ + approximate_inverse_->apply(alpha, b, beta, x); +} + + template void Isai::generate_inverse( std::shared_ptr input, bool skip_sorting, int power, @@ -361,6 +404,12 @@ std::unique_ptr Isai::conj_transpose() } +template +Isai::Isai(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + #define GKO_DECLARE_LOWER_ISAI(ValueType, IndexType) \ class Isai GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWER_ISAI); diff --git a/core/preconditioner/jacobi.cpp b/core/preconditioner/jacobi.cpp index 68a4d2d40d8..c5f8474f8d6 100644 --- a/core/preconditioner/jacobi.cpp +++ b/core/preconditioner/jacobi.cpp @@ -96,6 +96,19 @@ Jacobi::parse(const config::pnode& config, return params; } + +template +Jacobi::Jacobi(std::shared_ptr exec) + : LinOp(exec, dim<2>{}, precision_v), + num_blocks_{}, + blocks_(exec), + conditioning_(exec) +{ + parameters_.block_pointers.set_executor(exec); + parameters_.storage_optimization.block_wise.set_executor(exec); +} + + template Jacobi& Jacobi::operator=( const Jacobi& other) @@ -321,6 +334,67 @@ std::unique_ptr Jacobi::conj_transpose() const } +template +Jacobi::Jacobi(const Factory* factory, + std::shared_ptr system_matrix) + + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + parameters_{factory->get_parameters()}, + storage_scheme_{this->compute_storage_scheme( + parameters_.max_block_size, parameters_.max_block_stride)}, + num_blocks_{parameters_.block_pointers.get_size() - 1}, + blocks_(factory->get_executor(), + storage_scheme_.compute_storage_space( + parameters_.block_pointers.get_size() - 1)), + conditioning_(factory->get_executor()) +{ + parameters_.block_pointers.set_executor(this->get_executor()); + parameters_.storage_optimization.block_wise.set_executor( + this->get_executor()); + this->generate(system_matrix.get(), parameters_.skip_sorting); +} + + +template +block_interleaved_storage_scheme< + typename Jacobi::index_type> +Jacobi::compute_storage_scheme( + uint32 max_block_size, uint32 param_max_block_stride) + +{ + uint32 default_block_stride = 32; + // If the executor is hip, the warp size is 32 or 64 + if (auto hip_exec = std::dynamic_pointer_cast( + this->get_executor())) { + default_block_stride = hip_exec->get_warp_size(); + } + uint32 max_block_stride = default_block_stride; + if (param_max_block_stride != 0) { + // if parameter max_block_stride is not zero, set max_block_stride = + // param_max_block_stride + max_block_stride = param_max_block_stride; + if (this->get_executor() != this->get_executor()->get_master() && + max_block_stride != default_block_stride) { + // only support the default value on the gpu device + GKO_NOT_SUPPORTED(this); + } + } + if (parameters_.max_block_size > max_block_stride || + parameters_.max_block_size < 1) { + GKO_NOT_SUPPORTED(this); + } + const auto group_size = static_cast( + max_block_stride / get_superior_power(uint32{2}, max_block_size)); + const auto block_offset = max_block_size; + const auto block_stride = group_size * block_offset; + const auto group_offset = max_block_size * block_stride; + return {static_cast(block_offset), + static_cast(group_offset), + get_significant_bit(group_size)}; +} + + template void Jacobi::detect_blocks( const matrix::Csr* system_matrix) diff --git a/core/reorder/scaled_reordered.cpp b/core/reorder/scaled_reordered.cpp index 264122c0b8f..773618c052f 100644 --- a/core/reorder/scaled_reordered.cpp +++ b/core/reorder/scaled_reordered.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -15,6 +15,61 @@ namespace experimental { namespace reorder { +template +ScaledReordered::ScaledReordered( + std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v), + permutation_array_{exec} +{} + + +template +ScaledReordered::ScaledReordered( + const Factory* factory, std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), system_matrix->get_size(), + precision_v), + parameters_{factory->get_parameters()}, + permutation_array_{factory->get_executor()} +{ + // For now only support square matrices. + GKO_ASSERT_IS_SQUARE_MATRIX(system_matrix); + + auto exec = this->get_executor(); + + system_matrix_ = gko::clone(exec, system_matrix); + + // Scale the system matrix if scaling coefficients are provided + if (parameters_.row_scaling) { + GKO_ASSERT_EQUAL_DIMENSIONS(parameters_.row_scaling, system_matrix_); + row_scaling_ = parameters_.row_scaling; + row_scaling_->apply(system_matrix_, system_matrix_); + } + if (parameters_.col_scaling) { + GKO_ASSERT_EQUAL_DIMENSIONS(parameters_.col_scaling, system_matrix_); + col_scaling_ = parameters_.col_scaling; + col_scaling_->rapply(system_matrix_, system_matrix_); + } + + // If a reordering factory is provided, generate the reordering and + // permute the system matrix accordingly. + if (parameters_.reordering) { + auto reordering = parameters_.reordering->generate(system_matrix_); + permutation_array_ = reordering->get_permutation_array(); + system_matrix_ = as>(system_matrix_) + ->permute(&permutation_array_); + } + + // Generate the inner operator with the scaled and reordered system + // matrix. If none is provided, use the Identity. + if (parameters_.inner_operator) { + inner_operator_ = parameters_.inner_operator->generate(system_matrix_); + } else { + inner_operator_ = gko::matrix::Identity::create( + exec, this->get_size()[0]); + } +} + + template void ScaledReordered::apply_impl(const LinOp* b, LinOp* x) const diff --git a/core/solver/bicg.cpp b/core/solver/bicg.cpp index 4d13174a881..4d0912c5889 100644 --- a/core/solver/bicg.cpp +++ b/core/solver/bicg.cpp @@ -243,6 +243,23 @@ void Bicg::apply_dense_impl( } +template +Bicg::Bicg(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Bicg::Bicg(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template void Bicg::apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const diff --git a/core/solver/bicgstab.cpp b/core/solver/bicgstab.cpp index d4ff695dfc6..18c63ec396c 100644 --- a/core/solver/bicgstab.cpp +++ b/core/solver/bicgstab.cpp @@ -269,6 +269,23 @@ void Bicgstab::apply_impl(const LinOp* alpha, const LinOp* b, } +template +Bicgstab::Bicgstab(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Bicgstab::Bicgstab(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template int workspace_traits>::num_arrays(const Solver&) { diff --git a/core/solver/cb_gmres.cpp b/core/solver/cb_gmres.cpp index 2f2cd0bf823..570bf883bca 100644 --- a/core/solver/cb_gmres.cpp +++ b/core/solver/cb_gmres.cpp @@ -195,6 +195,23 @@ typename CbGmres::parameters_type CbGmres::parse( } +template +CbGmres::CbGmres(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +CbGmres::CbGmres(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template void CbGmres::apply_impl(const LinOp* b, LinOp* x) const { diff --git a/core/solver/cg.cpp b/core/solver/cg.cpp index 659e82fb1e0..743da453dc3 100644 --- a/core/solver/cg.cpp +++ b/core/solver/cg.cpp @@ -207,6 +207,25 @@ void Cg::apply_impl(const LinOp* alpha, const LinOp* b, } +template +Cg::Cg(std::shared_ptr exec) + + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Cg::Cg(const Factory* factory, + std::shared_ptr system_matrix) + + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template int workspace_traits>::num_arrays(const Solver&) { diff --git a/core/solver/cgs.cpp b/core/solver/cgs.cpp index 67680450842..01887f72d5a 100644 --- a/core/solver/cgs.cpp +++ b/core/solver/cgs.cpp @@ -239,6 +239,23 @@ int workspace_traits>::num_arrays(const Solver&) } +template +Cgs::Cgs(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Cgs::Cgs(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template int workspace_traits>::num_vectors(const Solver&) { diff --git a/core/solver/chebyshev.cpp b/core/solver/chebyshev.cpp index 348418a78d2..49fd676c58d 100644 --- a/core/solver/chebyshev.cpp +++ b/core/solver/chebyshev.cpp @@ -72,14 +72,15 @@ typename Chebyshev::parameters_type Chebyshev::parse( template Chebyshev::Chebyshev(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} template Chebyshev::Chebyshev(const Factory* factory, std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size())), + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), EnablePreconditionedIterativeSolver>{ std::move(system_matrix), factory->get_parameters()}, parameters_{factory->get_parameters()} diff --git a/core/solver/direct.cpp b/core/solver/direct.cpp index ede9ba46bbb..a01cf8672ce 100644 --- a/core/solver/direct.cpp +++ b/core/solver/direct.cpp @@ -52,7 +52,7 @@ std::unique_ptr Direct::conj_transpose() const template Direct::Direct(const Direct& other) - : LinOp{other.get_executor()} + : LinOp{other.get_executor(), dim<2>{}, precision_v} { *this = other; } @@ -60,7 +60,7 @@ Direct::Direct(const Direct& other) template Direct::Direct(Direct&& other) - : LinOp{other.get_executor()} + : LinOp{other.get_executor(), dim<2>{}, precision_v} { *this = std::move(other); } @@ -100,7 +100,7 @@ Direct& Direct::operator=( template Direct::Direct(std::shared_ptr exec) - : LinOp{exec} + : LinOp{exec, dim<2>{}, precision_v} {} @@ -124,7 +124,8 @@ generate_factorization( template Direct::Direct(const Factory* factory, std::shared_ptr system_matrix) - : LinOp{factory->get_executor(), system_matrix->get_size()}, + : LinOp{factory->get_executor(), system_matrix->get_size(), + precision_v}, gko::solver::EnableSolverBase< Direct, factorization::Factorization>{ generate_factorization( diff --git a/core/solver/fcg.cpp b/core/solver/fcg.cpp index 2b27131e499..801877f8f01 100644 --- a/core/solver/fcg.cpp +++ b/core/solver/fcg.cpp @@ -218,6 +218,23 @@ int workspace_traits>::num_arrays(const Solver&) } +template +Fcg::Fcg(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Fcg::Fcg(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template int workspace_traits>::num_vectors(const Solver&) { diff --git a/core/solver/gcr.cpp b/core/solver/gcr.cpp index 9c1bcba6cdf..15d7d9db8ca 100644 --- a/core/solver/gcr.cpp +++ b/core/solver/gcr.cpp @@ -298,6 +298,27 @@ void Gcr::apply_dense_impl(const VectorType* dense_b, } +template +Gcr::Gcr(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Gcr::Gcr(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{ + if (!parameters_.krylov_dim) { + parameters_.krylov_dim = gcr_default_krylov_dim; + } +} + + template void Gcr::apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const diff --git a/core/solver/gmres.cpp b/core/solver/gmres.cpp index a07602b7679..4bfb4f1e615 100644 --- a/core/solver/gmres.cpp +++ b/core/solver/gmres.cpp @@ -674,6 +674,27 @@ int workspace_traits>::num_arrays(const Solver&) } +template +Gmres::Gmres(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Gmres::Gmres(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{ + if (!parameters_.krylov_dim) { + parameters_.krylov_dim = gmres_default_krylov_dim; + } +} + + template int workspace_traits>::num_vectors(const Solver&) { diff --git a/core/solver/idr.cpp b/core/solver/idr.cpp index 612e449e5e2..51e959badf2 100644 --- a/core/solver/idr.cpp +++ b/core/solver/idr.cpp @@ -312,6 +312,23 @@ void Idr::iterate(const VectorType* dense_b, } +template +Idr::Idr(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Idr::Idr(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template void Idr::apply_impl(const LinOp* b, LinOp* x) const { diff --git a/core/solver/ir.cpp b/core/solver/ir.cpp index 7bfb42aaeaf..6a800837d4c 100644 --- a/core/solver/ir.cpp +++ b/core/solver/ir.cpp @@ -91,6 +91,37 @@ void Ir::set_relaxation_factor( } +template +Ir::Ir(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +Ir::Ir(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnableSolverBase{std::move(system_matrix)}, + EnableIterativeBase{ + stop::combine(factory->get_parameters().criteria)}, + parameters_{factory->get_parameters()} +{ + if (parameters_.generated_solver) { + this->set_solver(parameters_.generated_solver); + } else if (parameters_.solver) { + this->set_solver( + parameters_.solver->generate(this->get_system_matrix())); + } else { + this->set_solver(matrix::Identity::create( + this->get_executor(), this->get_size()[0])); + } + this->set_default_initial_guess(parameters_.default_initial_guess); + relaxation_factor_ = gko::initialize>( + {parameters_.relaxation_factor}, this->get_executor()); +} + + template Ir& Ir::operator=(const Ir& other) { diff --git a/core/solver/lower_trs.cpp b/core/solver/lower_trs.cpp index ccef06e8f84..bd83c539992 100644 --- a/core/solver/lower_trs.cpp +++ b/core/solver/lower_trs.cpp @@ -52,7 +52,7 @@ LowerTrs::parse( template LowerTrs::LowerTrs(const LowerTrs& other) - : LinOp(other.get_executor()) + : LinOp(other.get_executor(), dim<2>{}, precision_v) { *this = other; } @@ -60,7 +60,7 @@ LowerTrs::LowerTrs(const LowerTrs& other) template LowerTrs::LowerTrs(LowerTrs&& other) - : LinOp(other.get_executor()) + : LinOp(other.get_executor(), dim<2>{}, precision_v) { *this = std::move(other); } @@ -132,6 +132,26 @@ void LowerTrs::generate() } +template +LowerTrs::LowerTrs(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +LowerTrs::LowerTrs( + const Factory* factory, std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnableSolverBase, CsrMatrix>{ + copy_and_convert_to(factory->get_executor(), + system_matrix)}, + parameters_{factory->get_parameters()} +{ + this->generate(); +} + + static bool needs_transpose(std::shared_ptr exec) { bool result{}; diff --git a/core/solver/minres.cpp b/core/solver/minres.cpp index e8e992f8551..eb3e8adca6e 100644 --- a/core/solver/minres.cpp +++ b/core/solver/minres.cpp @@ -324,14 +324,15 @@ void Minres::apply_impl(const LinOp* alpha, const LinOp* b, template Minres::Minres(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} template Minres::Minres(const Factory* factory, std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size())), + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), EnablePreconditionedIterativeSolver{ std::move(system_matrix), factory->get_parameters()}, parameters_{factory->get_parameters()} diff --git a/core/solver/multigrid.cpp b/core/solver/multigrid.cpp index f1c5dbe114a..39f64ad2154 100644 --- a/core/solver/multigrid.cpp +++ b/core/solver/multigrid.cpp @@ -1073,7 +1073,8 @@ void Multigrid::create_state() const Multigrid::Multigrid(const Multigrid::Factory* factory, std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), transpose(system_matrix->get_size())), + : LinOp(factory->get_executor(), transpose(system_matrix->get_size()), + precision::any), EnableSolverBase{std::move(system_matrix)}, EnableIterativeBase{ stop::combine(factory->get_parameters().criteria)}, @@ -1107,7 +1108,9 @@ Multigrid::Multigrid(const Multigrid::Factory* factory, } -Multigrid::Multigrid(std::shared_ptr exec) : LinOp(exec) {} +Multigrid::Multigrid(std::shared_ptr exec) + : LinOp(exec, dim<2>{}, precision::any) +{} int workspace_traits::num_arrays(const Solver&) { return 1; } diff --git a/core/solver/pipe_cg.cpp b/core/solver/pipe_cg.cpp index 46d3787d188..34d3175f656 100644 --- a/core/solver/pipe_cg.cpp +++ b/core/solver/pipe_cg.cpp @@ -328,6 +328,23 @@ int workspace_traits>::num_arrays(const Solver&) } +template +PipeCg::PipeCg(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +PipeCg::PipeCg(const Factory* factory, + std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnablePreconditionedIterativeSolver>{ + std::move(system_matrix), factory->get_parameters()}, + parameters_{factory->get_parameters()} +{} + + template int workspace_traits>::num_vectors(const Solver&) { diff --git a/core/solver/upper_trs.cpp b/core/solver/upper_trs.cpp index 20d73c6fee6..f2fc56935e9 100644 --- a/core/solver/upper_trs.cpp +++ b/core/solver/upper_trs.cpp @@ -53,7 +53,7 @@ UpperTrs::parse( template UpperTrs::UpperTrs(const UpperTrs& other) - : LinOp(other.get_executor()) + : LinOp(other.get_executor(), dim<2>{}, precision_v) { *this = other; } @@ -61,7 +61,7 @@ UpperTrs::UpperTrs(const UpperTrs& other) template UpperTrs::UpperTrs(UpperTrs&& other) - : LinOp(other.get_executor()) + : LinOp(other.get_executor(), dim<2>{}, precision_v) { *this = std::move(other); } @@ -133,6 +133,26 @@ void UpperTrs::generate() } +template +UpperTrs::UpperTrs(std::shared_ptr exec) + : LinOp(std::move(exec), dim<2>{}, precision_v) +{} + + +template +UpperTrs::UpperTrs( + const Factory* factory, std::shared_ptr system_matrix) + : LinOp(factory->get_executor(), gko::transpose(system_matrix->get_size()), + precision_v), + EnableSolverBase, CsrMatrix>{ + copy_and_convert_to(factory->get_executor(), + system_matrix)}, + parameters_{factory->get_parameters()} +{ + this->generate(); +} + + static bool needs_transpose(std::shared_ptr exec) { bool result{}; diff --git a/core/test/base/CMakeLists.txt b/core/test/base/CMakeLists.txt index 09eb835e52e..0359923e912 100644 --- a/core/test/base/CMakeLists.txt +++ b/core/test/base/CMakeLists.txt @@ -32,5 +32,6 @@ ginkgo_create_test(sanitizers ADDITIONAL_LIBRARIES Threads::Threads) ginkgo_create_test(segmented_array LABELS distributed) ginkgo_create_test(segmented_range) ginkgo_create_test(types) +ginkgo_create_test(precision) ginkgo_create_test(utils) ginkgo_create_test(version EXECUTABLE_NAME version_test) # version collides with C++ stdlib header diff --git a/core/test/base/lin_op.cpp b/core/test/base/lin_op.cpp index 72ea7efc866..d7eef135c2c 100644 --- a/core/test/base/lin_op.cpp +++ b/core/test/base/lin_op.cpp @@ -2,6 +2,8 @@ // // SPDX-License-Identifier: BSD-3-Clause +#include "ginkgo/core/base/lin_op.hpp" + #include #include #include @@ -9,9 +11,10 @@ #include #include -#include #include +#include "core/test/utils.hpp" + namespace { @@ -75,8 +78,9 @@ class DummyLinOp : public gko::LinOp, public gko::EnableCreateMethod { public: DummyLinOp(std::shared_ptr exec, - gko::dim<2> size = gko::dim<2>{}) - : LinOp(exec, size) + gko::dim<2> size = gko::dim<2>{}, + gko::precision prec = gko::precision::none) + : LinOp(exec, size, prec) {} void access() const { last_access = this->get_executor(); } @@ -113,9 +117,33 @@ class DummyLinOp : public gko::LinOp, }; -class LinOp : public ::testing::Test { +TEST(LinOp, DefaultConstructs) +{ + auto exec = gko::ReferenceExecutor::create(); + auto op = std::make_shared(exec); + + EXPECT_EQ(op->get_executor(), exec); + EXPECT_EQ(op->get_size(), gko::dim<2>{}); + EXPECT_EQ(op->get_precision(), gko::precision::none); +} + + +TEST(LinOp, ConstructsWithArguments) +{ + auto exec = gko::ReferenceExecutor::create(); + auto size = gko::dim<2>{2, 3}; + auto prec = gko::precision::fp32; + auto op = std::make_shared(exec, size, prec); + + EXPECT_EQ(op->get_executor(), exec); + EXPECT_EQ(op->get_size(), size); + EXPECT_EQ(op->get_precision(), prec); +} + + +class LinOpApply : public ::testing::Test { protected: - LinOp() + LinOpApply() : ref{gko::ReferenceExecutor::create()}, ref2{gko::ReferenceExecutor::create()}, op{DummyLinOp::create(ref2, gko::dim<2>{3, 5})}, @@ -139,7 +167,7 @@ class LinOp : public ::testing::Test { }; -TEST_F(LinOp, CallsApplyImpl) +TEST_F(LinOpApply, CallsImpl) { op->apply(b, x); @@ -147,7 +175,7 @@ TEST_F(LinOp, CallsApplyImpl) } -TEST_F(LinOp, CallsExtendedApplyImpl) +TEST_F(LinOpApply, CallsExtendedImpl) { op->apply(alpha, b, beta, x); @@ -155,7 +183,7 @@ TEST_F(LinOp, CallsExtendedApplyImpl) } -TEST_F(LinOp, ApplyFailsOnWrongBSize) +TEST_F(LinOpApply, FailsOnWrongBSize) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{3, 4}); @@ -163,7 +191,7 @@ TEST_F(LinOp, ApplyFailsOnWrongBSize) } -TEST_F(LinOp, ApplyFailsOnWrongSolutionRows) +TEST_F(LinOpApply, FailsOnWrongSolutionRows) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{5, 4}); @@ -171,7 +199,7 @@ TEST_F(LinOp, ApplyFailsOnWrongSolutionRows) } -TEST_F(LinOp, ApplyFailsOnWrongSolutionColumns) +TEST_F(LinOpApply, FailsOnWrongSolutionColumns) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{3, 5}); @@ -179,7 +207,7 @@ TEST_F(LinOp, ApplyFailsOnWrongSolutionColumns) } -TEST_F(LinOp, ExtendedApplyFailsOnWrongBSize) +TEST_F(LinOpApply, ExtendedFailsOnWrongBSize) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{3, 4}); @@ -187,7 +215,7 @@ TEST_F(LinOp, ExtendedApplyFailsOnWrongBSize) } -TEST_F(LinOp, ExtendedApplyFailsOnWrongSolutionRows) +TEST_F(LinOpApply, ExtendedFailsOnWrongSolutionRows) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{5, 4}); @@ -195,7 +223,7 @@ TEST_F(LinOp, ExtendedApplyFailsOnWrongSolutionRows) } -TEST_F(LinOp, ExtendedApplyFailsOnWrongSolutionColumns) +TEST_F(LinOpApply, ExtendedFailsOnWrongSolutionColumns) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{3, 5}); @@ -203,7 +231,7 @@ TEST_F(LinOp, ExtendedApplyFailsOnWrongSolutionColumns) } -TEST_F(LinOp, ExtendedApplyFailsOnWrongAlphaDimension) +TEST_F(LinOpApply, ExtendedFailsOnWrongAlphaDimension) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{2, 5}); @@ -211,7 +239,7 @@ TEST_F(LinOp, ExtendedApplyFailsOnWrongAlphaDimension) } -TEST_F(LinOp, ExtendedApplyFailsOnWrongBetaDimension) +TEST_F(LinOpApply, ExtendedFailsOnWrongBetaDimension) { auto wrong = DummyLinOp::create(ref, gko::dim<2>{2, 5}); @@ -220,7 +248,7 @@ TEST_F(LinOp, ExtendedApplyFailsOnWrongBetaDimension) // For tests between different memory, check cuda/test/base/lin_op.cu -TEST_F(LinOp, ApplyDoesNotCopyBetweenSameMemory) +TEST_F(LinOpApply, DoesNotCopyBetweenSameMemory) { op->apply(b, x); @@ -229,7 +257,7 @@ TEST_F(LinOp, ApplyDoesNotCopyBetweenSameMemory) } -TEST_F(LinOp, ApplyNoCopyBackBetweenSameMemory) +TEST_F(LinOpApply, NoCopyBackBetweenSameMemory) { op->apply(b, x); @@ -238,7 +266,7 @@ TEST_F(LinOp, ApplyNoCopyBackBetweenSameMemory) } -TEST_F(LinOp, ExtendedApplyDoesNotCopyBetweenSameMemory) +TEST_F(LinOpApply, ExtendedDoesNotCopyBetweenSameMemory) { op->apply(alpha, b, beta, x); @@ -249,7 +277,7 @@ TEST_F(LinOp, ExtendedApplyDoesNotCopyBetweenSameMemory) } -TEST_F(LinOp, ExtendedApplyNoCopyBackBetweenSameMemory) +TEST_F(LinOpApply, ExtendedNoCopyBackBetweenSameMemory) { op->apply(alpha, b, beta, x); @@ -260,13 +288,13 @@ TEST_F(LinOp, ExtendedApplyNoCopyBackBetweenSameMemory) } -TEST_F(LinOp, ApplyUsesInitialGuessReturnsFalse) +TEST_F(LinOpApply, UsesInitialGuessReturnsFalse) { ASSERT_FALSE(op->apply_uses_initial_guess()); } -TEST_F(LinOp, ApplyIsLogged) +TEST_F(LinOpApply, IsLogged) { auto before_logger = *logger; @@ -279,7 +307,7 @@ TEST_F(LinOp, ApplyIsLogged) } -TEST_F(LinOp, AdvancedApplyIsLogged) +TEST_F(LinOpApply, AdvancedIsLogged) { auto before_logger = *logger; @@ -430,7 +458,7 @@ class DummyLinOpWithType DummyLinOpWithType(std::shared_ptr exec, gko::dim<2> size, Type value) - : gko::LinOp(exec, size), value_(value) + : gko::LinOp(exec, size, gko::precision_v), value_(value) {} void compute_absolute_inplace() override { value_ = gko::abs(value_); } @@ -455,6 +483,33 @@ class DummyLinOpWithType }; +template +class LinOpPrecision : public ::testing::Test { +protected: + using value_type = ValueType; + using dummy_type = DummyLinOpWithType; + + LinOpPrecision() + : ref{gko::ReferenceExecutor::create()}, + op{dummy_type::create(ref, gko::dim<2>{1, 1}, value_type{})} + {} + + std::shared_ptr ref; + std::shared_ptr op; +}; + +TYPED_TEST_SUITE(LinOpPrecision, gko::test::ValueTypes, TypenameNameGenerator); + + +TYPED_TEST(LinOpPrecision, HasCorrectPrecision) +{ + using value_type = typename TestFixture::value_type; + auto expected_precision = gko::precision_v; + + EXPECT_EQ(expected_precision, this->op->get_precision()); +} + + class EnableAbsoluteComputation : public ::testing::Test { protected: using dummy_type = DummyLinOpWithType>; diff --git a/core/test/base/precision.cpp b/core/test/base/precision.cpp new file mode 100644 index 00000000000..a7bb64152fe --- /dev/null +++ b/core/test/base/precision.cpp @@ -0,0 +1,181 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "ginkgo/core/base/precision.hpp" + +#include + +#include + +namespace { +namespace gko_test { + + +auto precisions = { + gko::precision::none, gko::precision::any, + gko::precision::fp32, gko::precision::complex_fp32, + gko::precision::fp64, gko::precision::complex_fp64, +#if GINKGO_ENABLE_BFLOAT16 + gko::precision::bf16, gko::precision::complex_bf16, +#endif +#if GINKGO_ENABLE_HALF + gko::precision::fp16, gko::precision::complex_fp16, +#endif +}; + + +TEST(Precision, EnumOpsEqual) +{ + auto any_p = gko::precision::any; + + for (auto p1 : precisions) { + auto same_p = p1; + + EXPECT_EQ(p1, same_p); + EXPECT_EQ(same_p, p1); + EXPECT_EQ(p1, any_p); + EXPECT_EQ(any_p, p1); + } +} + + +TEST(Precision, EnumOpsNotEqual) +{ + for (auto p1 : precisions) { + for (auto p2 : precisions) { + if (p1 == p2) { + continue; + } + + EXPECT_NE(p1, p2); + EXPECT_NE(p2, p1); + } + } +} + + +TEST(Precision, TypeToPrecision) +{ + EXPECT_EQ(gko::precision_v, gko::precision::fp32); + EXPECT_EQ(gko::precision_v>, + gko::precision::complex_fp32); + EXPECT_EQ(gko::precision_v, gko::precision::fp64); + EXPECT_EQ(gko::precision_v>, + gko::precision::complex_fp64); +#if GINKGO_ENABLE_BFLOAT16 + EXPECT_EQ(gko::precision_v, gko::precision::bf16); + EXPECT_EQ(gko::precision_v>, + gko::precision::complex_bf16); +#endif +#if GINKGO_ENABLE_HALF + EXPECT_EQ(gko::precision_v, gko::precision::fp16); + EXPECT_EQ(gko::precision_v>, + gko::precision::complex_fp16); +#endif +} + + +template +void test_is_complex() +{ + EXPECT_EQ(gko::is_complex(), gko::is_complex(gko::precision_v)); +} + +TEST(Precision, IsComplex) +{ + test_is_complex(); + test_is_complex(); + test_is_complex>(); + test_is_complex>(); +#if GINKGO_ENABLE_HALF + test_is_complex(); + test_is_complex>(); +#endif +#if GINKGO_ENABLE_BFLOAT16 + test_is_complex(); + test_is_complex>(); +#endif +} + + +template +void test_is_real() +{ + EXPECT_EQ(!gko::is_complex(), gko::is_real(gko::precision_v)); +} + +TEST(Precision, IsReal) +{ + test_is_real(); + test_is_real(); + test_is_real>(); + test_is_real>(); +#if GINKGO_ENABLE_HALF + test_is_real(); + test_is_real>(); +#endif +#if GINKGO_ENABLE_BFLOAT16 + test_is_real(); + test_is_real>(); +#endif +} + +TEST(Precision, AsReal) +{ + EXPECT_EQ(gko::as_real(gko::precision::fp32), gko::precision::fp32); + EXPECT_EQ(gko::as_real(gko::precision::fp64), gko::precision::fp64); + EXPECT_EQ(gko::as_real(gko::precision::complex_fp32), gko::precision::fp32); + EXPECT_EQ(gko::as_real(gko::precision::complex_fp64), gko::precision::fp64); + EXPECT_EQ(gko::as_real(gko::precision::complex_fp16), gko::precision::fp16); + EXPECT_EQ(gko::as_real(gko::precision::any), gko::precision::any); + EXPECT_THROW(gko::as_real(gko::precision::none), gko::InvalidStateError); +} + +TEST(Precision, AsComplex) +{ + EXPECT_EQ(gko::as_complex(gko::precision::fp32), + gko::precision::complex_fp32); + EXPECT_EQ(gko::as_complex(gko::precision::fp64), + gko::precision::complex_fp64); + EXPECT_EQ(gko::as_complex(gko::precision::complex_fp32), + gko::precision::complex_fp32); + EXPECT_EQ(gko::as_complex(gko::precision::complex_fp64), + gko::precision::complex_fp64); + EXPECT_EQ(gko::as_complex(gko::precision::fp16), + gko::precision::complex_fp16); + EXPECT_EQ(gko::as_real(gko::precision::any), gko::precision::any); + EXPECT_THROW(gko::as_complex(gko::precision::none), gko::InvalidStateError); +} + +TEST(Precision, PrecisionToVariant) +{ + EXPECT_TRUE(std::holds_alternative( + gko::precision_to_variant(gko::precision::fp32))); + EXPECT_TRUE(std::holds_alternative( + gko::precision_to_variant(gko::precision::fp64))); + EXPECT_TRUE(std::holds_alternative>( + gko::precision_to_variant(gko::precision::complex_fp32))); + EXPECT_TRUE(std::holds_alternative>( + gko::precision_to_variant(gko::precision::complex_fp64))); +#if GINKGO_ENABLE_HALF + EXPECT_TRUE(std::holds_alternative( + gko::precision_to_variant(gko::precision::fp16))); + EXPECT_TRUE(std::holds_alternative>( + gko::precision_to_variant(gko::precision::complex_fp16))); +#endif +#if GINKGO_ENABLE_BFLOAT16 + EXPECT_TRUE(std::holds_alternative( + gko::precision_to_variant(gko::precision::bf16))); + EXPECT_TRUE(std::holds_alternative>( + gko::precision_to_variant(gko::precision::complex_bf16))); +#endif + EXPECT_THROW(gko::precision_to_variant(gko::precision::none), + gko::InvalidStateError); + EXPECT_THROW(gko::precision_to_variant(gko::precision::any), + gko::InvalidStateError); +} + + +} // namespace gko_test +} // namespace diff --git a/core/test/matrix/dense.cpp b/core/test/matrix/dense.cpp index 8e04d43e0a4..706c84405ce 100644 --- a/core/test/matrix/dense.cpp +++ b/core/test/matrix/dense.cpp @@ -40,7 +40,7 @@ class Dense : public ::testing::Test { static void assert_empty(gko::ptr_param> m) { - ASSERT_EQ(m->get_precision(), gko::type_to_precision); + ASSERT_EQ(m->get_precision(), gko::precision_v); ASSERT_EQ(m->get_size(), gko::dim<2>(0, 0)); ASSERT_EQ(m->get_num_stored_elements(), 0); } diff --git a/extensions/cuda/solver/cudss.cpp b/extensions/cuda/solver/cudss.cpp index e287f7a661c..171118d31d8 100644 --- a/extensions/cuda/solver/cudss.cpp +++ b/extensions/cuda/solver/cudss.cpp @@ -94,13 +94,13 @@ struct Cudss::state { template Cudss::Cudss(std::shared_ptr exec) - : LinOp{exec} + : LinOp{exec, dim<2>{}, precision_v} {} template Cudss::Cudss(const Cudss& other) - : LinOp{other.get_executor()} + : Cudss{other.get_executor()} { *this = other; } @@ -108,7 +108,7 @@ Cudss::Cudss(const Cudss& other) template Cudss::Cudss(Cudss&& other) noexcept - : LinOp{other.get_executor()} + : Cudss{other.get_executor()} { *this = std::move(other); } @@ -143,7 +143,8 @@ Cudss& Cudss::operator=( template Cudss::Cudss(const Factory* factory, std::shared_ptr system_matrix) - : LinOp{factory->get_executor(), system_matrix->get_size()} + : LinOp{factory->get_executor(), system_matrix->get_size(), + precision_v} { const auto exec = this->get_executor(); auto cuda_exec = std::dynamic_pointer_cast(exec); diff --git a/include/ginkgo/core/base/combination.hpp b/include/ginkgo/core/base/combination.hpp index aeb32045dd9..32d22487845 100644 --- a/include/ginkgo/core/base/combination.hpp +++ b/include/ginkgo/core/base/combination.hpp @@ -117,7 +117,9 @@ class Combination : public LinOp, * * @param exec Executor associated to the linear combination */ - explicit Combination(std::shared_ptr exec) : LinOp(exec) {} + explicit Combination(std::shared_ptr exec) + : LinOp(exec, dim<2>{}, precision_v) + {} /** * Creates a linear combination of operators using the specified list of @@ -143,7 +145,7 @@ class Combination : public LinOp, CoefficientIterator coefficient_end, OperatorIterator operator_begin, OperatorIterator operator_end) - : LinOp([&] { + : Combination([&] { if (operator_begin == operator_end) { throw OutOfBoundsError(__FILE__, __LINE__, 1, 0); } diff --git a/include/ginkgo/core/base/composition.hpp b/include/ginkgo/core/base/composition.hpp index 156761c0d87..93b46cf8eee 100644 --- a/include/ginkgo/core/base/composition.hpp +++ b/include/ginkgo/core/base/composition.hpp @@ -113,7 +113,7 @@ class Composition : public LinOp, * @param exec Executor associated to the composition */ explicit Composition(std::shared_ptr exec) - : LinOp(exec), storage_{exec} + : LinOp(exec, dim<2>{}, precision_v), storage_{exec} {} /** @@ -129,13 +129,12 @@ class Composition : public LinOp, typename = std::void_t< typename std::iterator_traits::iterator_category>> explicit Composition(Iterator begin, Iterator end) - : LinOp([&] { + : Composition([&] { if (begin == end) { throw OutOfBoundsError(__FILE__, __LINE__, 1, 0); } return (*begin)->get_executor(); - }()), - storage_{this->get_executor()} + }()) { for (auto it = begin; it != end; ++it) { add_operators(*it); diff --git a/include/ginkgo/core/base/exception.hpp b/include/ginkgo/core/base/exception.hpp index 00661052027..844fb948e9d 100644 --- a/include/ginkgo/core/base/exception.hpp +++ b/include/ginkgo/core/base/exception.hpp @@ -11,6 +11,8 @@ #include +#include "precision.hpp" + namespace gko { @@ -673,6 +675,30 @@ class InvalidStateError : public Error { }; +/** + * Exception thrown if an object with mismatched precision is used + */ +class PrecisionError : public Error { +public: + /** + * + * @param file The name of the offending source file + * @param line The source code line number where the error occurred + * @param func The function name where the error occurred + * @param expected The string value of the expected precision + * @param op The offending operator + * @param p The string value of the precision of the operator + */ + PrecisionError(const std::string& file, int line, const std::string& func, + precision expected, const std::string& op, precision p) + : Error(file, line, + func + ": expected precision " + to_string(expected) + + " but object " + op + " has precision " + to_string(p) + + ".") + {} +}; + + } // namespace gko diff --git a/include/ginkgo/core/base/exception_helpers.hpp b/include/ginkgo/core/base/exception_helpers.hpp index f0104ba1a7c..678f298ea22 100644 --- a/include/ginkgo/core/base/exception_helpers.hpp +++ b/include/ginkgo/core/base/exception_helpers.hpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2026 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -12,6 +12,7 @@ #include #include #include +#include #include @@ -119,6 +120,15 @@ inline dim<2> get_size(const T& op) inline dim<2> get_size(const dim<2>& size) { return size; } +template +precision get_precision(T&& op) +{ + return op->get_precision(); +} + +inline precision get_precision(precision p) { return p; } + + template inline batch_dim<2> get_batch_size(const T& op) { @@ -286,6 +296,21 @@ inline size_type get_num_batch_items(const T& obj) } +/** + * Asserts that _op has the expected precision _expected. + * + * @throw PrecisionError if _op has a precision other than _expected. + */ +#define GKO_ASSERT_PRECISION(_expected, _op) \ + if (::gko::detail::get_precision(_expected) != \ + ::gko::detail::get_precision(_op)) { \ + throw ::gko::PrecisionError(__FILE__, __LINE__, __func__, \ + ::gko::detail::get_precision(_expected), \ + #_op, ::gko::detail::get_precision(_op)); \ + } \ + static_assert(true, "Require ;") + + /** * Asserts that _op1 and _op2 have equal number of items in the batch * diff --git a/include/ginkgo/core/base/lin_op.hpp b/include/ginkgo/core/base/lin_op.hpp index 5e87a61a9a3..232a930c200 100644 --- a/include/ginkgo/core/base/lin_op.hpp +++ b/include/ginkgo/core/base/lin_op.hpp @@ -22,6 +22,8 @@ #include #include +#include "precision.hpp" + namespace gko { namespace matrix { @@ -177,6 +179,8 @@ class LinOp : public PolymorphicObject { */ virtual bool apply_uses_initial_guess() const { return false; } + [[nodiscard]] precision get_precision() const noexcept; + /** Copy-assigns a LinOp. Preserves the executor and copies the size. */ LinOp& operator=(const LinOp&) = default; @@ -202,10 +206,7 @@ class LinOp : public PolymorphicObject { * Move-constructs a LinOp. Inherits executor and size from the input, * which will have size 0x0 and unchanged executor afterwards. */ - LinOp(LinOp&& other) - : PolymorphicObject(std::move(other)), - size_{std::exchange(other.size_, dim<2>{})} - {} + LinOp(LinOp&& other); protected: /** @@ -215,9 +216,8 @@ class LinOp : public PolymorphicObject { * @param size the size of the operator */ explicit LinOp(std::shared_ptr exec, - const dim<2>& size = dim<2>{}) - : PolymorphicObject(exec), size_{size} - {} + const dim<2>& size = dim<2>{}, + precision p = precision::none); /** * Sets the size of the operator. @@ -226,6 +226,8 @@ class LinOp : public PolymorphicObject { */ void set_size(const dim<2>& value) noexcept { size_ = value; } + void set_precision(precision p) noexcept; + /** * Implementers of LinOp should override this function instead * of apply(const LinOp *, LinOp *). @@ -283,6 +285,7 @@ class LinOp : public PolymorphicObject { private: dim<2> size_{}; + precision value_t_{precision::none}; }; @@ -854,7 +857,8 @@ using EnableDefaultLinOpFactory = * // constructor needed by the factory * explicit MyLinOp(const Factory *factory, * std::shared_ptr matrix) - * : LinOp(factory->get_executor()), matrix->get_size()), + * : LinOp(factory->get_executor()), matrix->get_size(), + * precision::float64), * // store factory's parameters locally * my_parameters_{factory->get_parameters()}, * { diff --git a/include/ginkgo/core/base/precision.hpp b/include/ginkgo/core/base/precision.hpp new file mode 100644 index 00000000000..c4c3ab36877 --- /dev/null +++ b/include/ginkgo/core/base/precision.hpp @@ -0,0 +1,171 @@ +// SPDX-FileCopyrightText: 2026 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#pragma once + +#include +#include + +#include +#include + +namespace gko { + + +/** + * A enum to specify the precision of stored data. + */ +enum struct precision { + none, //!< no precision information is available, incompatible with all + //!< other precisions + any, //!< compatible with all other precisions, including none + fp32, + complex_fp32, + fp64, + complex_fp64, +#if GINKGO_ENABLE_HALF + fp16, + complex_fp16, +#endif +#if GINKGO_ENABLE_BFLOAT16 + bf16, + complex_bf16, +#endif +}; + + +// Returns the string representation of a precision +std::string to_string(precision p); + + +// Equality comparison operator. +// The precision::any is equal to any other precision, including none. +constexpr bool operator==(precision a, precision b) +{ + auto int_a = static_cast(a); + auto int_b = static_cast(b); + if (int_a == static_cast(precision::any) || + int_b == static_cast(precision::any)) { + return true; + } + return int_a == int_b; +} + + +constexpr bool operator!=(precision a, precision b) { return !(a == b); } + + +/** + * Map from compile time type to runtime precision. + * + * @tparam T Value type to map to a precision + */ +template +inline precision precision_v; + +template <> +inline constexpr precision precision_v = precision::fp32; +template <> +inline constexpr precision precision_v> = + precision::complex_fp32; +template <> +inline constexpr precision precision_v = precision::fp64; +template <> +inline constexpr precision precision_v> = + precision::complex_fp64; +#if GINKGO_ENABLE_HALF +template <> +inline constexpr precision precision_v = precision::fp16; +template <> +inline constexpr precision precision_v> = + precision::complex_fp16; +#endif +#if GINKGO_ENABLE_BFLOAT16 +template <> +inline constexpr precision precision_v = precision::bf16; +template <> +inline constexpr precision precision_v> = + precision::complex_bf16; +#endif + + +// True if the precision is complex or any +constexpr bool is_complex(precision p) +{ + return +#if GINKGO_ENABLE_HALF + p == precision::complex_fp16 || +#endif +#if GINKGO_ENABLE_BFLOAT16 + p == precision::complex_bf16 || +#endif + p == precision::complex_fp32 || p == precision::complex_fp64; +} + + +// True if the precision is real or any +constexpr bool is_real(precision p) +{ + return +#if GINKGO_ENABLE_HALF + p == precision::fp16 || +#endif +#if GINKGO_ENABLE_BFLOAT16 + p == precision::bf16 || +#endif + p == precision::fp32 || p == precision::fp64; +} + + +/** + * Maps a precision to its corresponding real precision. + * + * For example as_real(complex_fp32) == fp32. + * For a real precision or any this is the identity. + * + * @throws InvalidStateError if the precision is none + * + * @param p precision to map to a real precision + * @return The real precision corresponding to p + */ +precision as_real(precision p); + + +/** + * Maps a precision to its corresponding complex precision. + * + * For example as_complex(fp32) == complex_fp32. + * For a complex precision or any this is the identity. + * + * @throws InvalidStateError if the precision is none + * + * @param p precision to map to a complex precision + * @return The complex precision corresponding to p + */ +precision as_complex(precision p); + + +/** + * Create a variant from a precision. + * + * This allows to map the runtime precision back to a compile time precision + * when needed. + * For example, for a precision fp32, the result variant will hold float as + * its alternatives. + * + * @param p The precision to map to the variant + * @return A variant which value corresponds to the type matching the precision. + */ +std::variant< +#if GINKGO_ENABLE_HALF + half, std::complex, +#endif +#if GINKGO_ENABLE_BFLOAT16 + bfloat16, std::complex, +#endif + float, std::complex, double, std::complex> +precision_to_variant(precision p); + + +} // namespace gko diff --git a/include/ginkgo/core/distributed/preconditioner/schwarz.hpp b/include/ginkgo/core/distributed/preconditioner/schwarz.hpp index 232a3e9843b..d2b447460d6 100644 --- a/include/ginkgo/core/distributed/preconditioner/schwarz.hpp +++ b/include/ginkgo/core/distributed/preconditioner/schwarz.hpp @@ -168,7 +168,7 @@ class Schwarz : public LinOp { * @param exec the executor this object is assigned to */ explicit Schwarz(std::shared_ptr exec) - : LinOp(std::move(exec)) + : LinOp(std::move(exec), dim<2>{}, precision_v) {} /** @@ -181,7 +181,8 @@ class Schwarz : public LinOp { explicit Schwarz(const Factory* factory, std::shared_ptr system_matrix) : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), + gko::transpose(system_matrix->get_size()), + precision_v), parameters_{factory->get_parameters()}, system_matrix_{system_matrix} { diff --git a/include/ginkgo/core/preconditioner/isai.hpp b/include/ginkgo/core/preconditioner/isai.hpp index d02c566dc18..54522806450 100644 --- a/include/ginkgo/core/preconditioner/isai.hpp +++ b/include/ginkgo/core/preconditioner/isai.hpp @@ -215,8 +215,7 @@ class Isai : public LinOp, public Transposable { std::unique_ptr conj_transpose() const override; protected: - explicit Isai(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Isai(std::shared_ptr exec); /** * Creates an Isai preconditioner from a matrix using an Isai::Factory. @@ -225,34 +224,12 @@ class Isai : public LinOp, public Transposable { * @param system_matrix the matrix for which an ISAI is to be computed */ explicit Isai(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), system_matrix->get_size()), - parameters_{factory->get_parameters()} - { - const auto skip_sorting = parameters_.skip_sorting; - const auto power = parameters_.sparsity_power; - const auto excess_limit = parameters_.excess_limit; - generate_inverse(system_matrix, skip_sorting, power, excess_limit, - static_cast>( - parameters_.excess_solver_reduction)); - if (IsaiType == isai_type::spd) { - auto inv = share(as(approximate_inverse_)); - auto inv_transp = share(inv->conj_transpose()); - approximate_inverse_ = - Composition::create(inv_transp, inv); - } - } + std::shared_ptr system_matrix); - void apply_impl(const LinOp* b, LinOp* x) const override - { - approximate_inverse_->apply(b, x); - } + void apply_impl(const LinOp* b, LinOp* x) const override; void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, - LinOp* x) const override - { - approximate_inverse_->apply(alpha, b, beta, x); - } + LinOp* x) const override; private: /** diff --git a/include/ginkgo/core/preconditioner/jacobi.hpp b/include/ginkgo/core/preconditioner/jacobi.hpp index 0b24558abd3..34ac73508c9 100644 --- a/include/ginkgo/core/preconditioner/jacobi.hpp +++ b/include/ginkgo/core/preconditioner/jacobi.hpp @@ -540,12 +540,7 @@ class Jacobi : public LinOp, * * @param exec the executor this object is assigned to */ - explicit Jacobi(std::shared_ptr exec) - : LinOp(exec), num_blocks_{}, blocks_(exec), conditioning_(exec) - { - parameters_.block_pointers.set_executor(exec); - parameters_.storage_optimization.block_wise.set_executor(exec); - } + explicit Jacobi(std::shared_ptr exec); /** * Creates a Jacobi preconditioner from a matrix using a Jacobi::Factory. @@ -555,23 +550,7 @@ class Jacobi : public LinOp, * from */ explicit Jacobi(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - parameters_{factory->get_parameters()}, - storage_scheme_{this->compute_storage_scheme( - parameters_.max_block_size, parameters_.max_block_stride)}, - num_blocks_{parameters_.block_pointers.get_size() - 1}, - blocks_(factory->get_executor(), - storage_scheme_.compute_storage_space( - parameters_.block_pointers.get_size() - 1)), - conditioning_(factory->get_executor()) - { - parameters_.block_pointers.set_executor(this->get_executor()); - parameters_.storage_optimization.block_wise.set_executor( - this->get_executor()); - this->generate(system_matrix.get(), parameters_.skip_sorting); - } + std::shared_ptr system_matrix); /** * Computes the storage scheme suitable for storing blocks of a given @@ -582,38 +561,7 @@ class Jacobi : public LinOp, * @return a suitable storage scheme */ block_interleaved_storage_scheme compute_storage_scheme( - uint32 max_block_size, uint32 param_max_block_stride) - { - uint32 default_block_stride = 32; - // If the executor is hip, the warp size is 32 or 64 - if (auto hip_exec = std::dynamic_pointer_cast( - this->get_executor())) { - default_block_stride = hip_exec->get_warp_size(); - } - uint32 max_block_stride = default_block_stride; - if (param_max_block_stride != 0) { - // if parameter max_block_stride is not zero, set max_block_stride = - // param_max_block_stride - max_block_stride = param_max_block_stride; - if (this->get_executor() != this->get_executor()->get_master() && - max_block_stride != default_block_stride) { - // only support the default value on the gpu device - GKO_NOT_SUPPORTED(this); - } - } - if (parameters_.max_block_size > max_block_stride || - parameters_.max_block_size < 1) { - GKO_NOT_SUPPORTED(this); - } - const auto group_size = static_cast( - max_block_stride / get_superior_power(uint32{2}, max_block_size)); - const auto block_offset = max_block_size; - const auto block_stride = group_size * block_offset; - const auto group_offset = max_block_size * block_stride; - return {static_cast(block_offset), - static_cast(group_offset), - get_significant_bit(group_size)}; - } + uint32 max_block_size, uint32 param_max_block_stride); /** * Generates the preconditioner. diff --git a/include/ginkgo/core/reorder/scaled_reordered.hpp b/include/ginkgo/core/reorder/scaled_reordered.hpp index dffc62653a1..8bc7831fe8b 100644 --- a/include/ginkgo/core/reorder/scaled_reordered.hpp +++ b/include/ginkgo/core/reorder/scaled_reordered.hpp @@ -97,56 +97,10 @@ class ScaledReordered : public LinOp { /** * Creates an empty scaled reordered operator (0x0 operator). */ - explicit ScaledReordered(std::shared_ptr exec) - : LinOp(std::move(exec)), permutation_array_{exec} - {} + explicit ScaledReordered(std::shared_ptr exec); explicit ScaledReordered(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), system_matrix->get_size()), - parameters_{factory->get_parameters()}, - permutation_array_{factory->get_executor()} - { - // For now only support square matrices. - GKO_ASSERT_IS_SQUARE_MATRIX(system_matrix); - - auto exec = this->get_executor(); - - system_matrix_ = gko::clone(exec, system_matrix); - - // Scale the system matrix if scaling coefficients are provided - if (parameters_.row_scaling) { - GKO_ASSERT_EQUAL_DIMENSIONS(parameters_.row_scaling, - system_matrix_); - row_scaling_ = parameters_.row_scaling; - row_scaling_->apply(system_matrix_, system_matrix_); - } - if (parameters_.col_scaling) { - GKO_ASSERT_EQUAL_DIMENSIONS(parameters_.col_scaling, - system_matrix_); - col_scaling_ = parameters_.col_scaling; - col_scaling_->rapply(system_matrix_, system_matrix_); - } - - // If a reordering factory is provided, generate the reordering and - // permute the system matrix accordingly. - if (parameters_.reordering) { - auto reordering = parameters_.reordering->generate(system_matrix_); - permutation_array_ = reordering->get_permutation_array(); - system_matrix_ = as>(system_matrix_) - ->permute(&permutation_array_); - } - - // Generate the inner operator with the scaled and reordered system - // matrix. If none is provided, use the Identity. - if (parameters_.inner_operator) { - inner_operator_ = - parameters_.inner_operator->generate(system_matrix_); - } else { - inner_operator_ = gko::matrix::Identity::create( - exec, this->get_size()[0]); - } - } + std::shared_ptr system_matrix); void apply_impl(const LinOp* b, LinOp* x) const override; diff --git a/include/ginkgo/core/solver/bicg.hpp b/include/ginkgo/core/solver/bicg.hpp index 04d55cc9749..68f044af555 100644 --- a/include/ginkgo/core/solver/bicg.hpp +++ b/include/ginkgo/core/solver/bicg.hpp @@ -107,17 +107,10 @@ class Bicg void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Bicg(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Bicg(std::shared_ptr exec); explicit Bicg(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/bicgstab.hpp b/include/ginkgo/core/solver/bicgstab.hpp index 47d3747b71f..0392d4484cb 100644 --- a/include/ginkgo/core/solver/bicgstab.hpp +++ b/include/ginkgo/core/solver/bicgstab.hpp @@ -104,18 +104,10 @@ class Bicgstab void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Bicgstab(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit Bicgstab(std::shared_ptr exec); explicit Bicgstab(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/cb_gmres.hpp b/include/ginkgo/core/solver/cb_gmres.hpp index 54ed00fbc37..ffaf0bc59b5 100644 --- a/include/ginkgo/core/solver/cb_gmres.hpp +++ b/include/ginkgo/core/solver/cb_gmres.hpp @@ -172,17 +172,10 @@ class CbGmres : public LinOp, void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit CbGmres(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit CbGmres(std::shared_ptr exec); explicit CbGmres(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/cg.hpp b/include/ginkgo/core/solver/cg.hpp index ae5130da591..ecee57bcd5f 100644 --- a/include/ginkgo/core/solver/cg.hpp +++ b/include/ginkgo/core/solver/cg.hpp @@ -100,17 +100,10 @@ class Cg : public LinOp, void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Cg(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Cg(std::shared_ptr exec); explicit Cg(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/cgs.hpp b/include/ginkgo/core/solver/cgs.hpp index 3a11b8cf23e..492b40469a9 100644 --- a/include/ginkgo/core/solver/cgs.hpp +++ b/include/ginkgo/core/solver/cgs.hpp @@ -97,17 +97,10 @@ class Cgs void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Cgs(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Cgs(std::shared_ptr exec); explicit Cgs(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/fcg.hpp b/include/ginkgo/core/solver/fcg.hpp index a713bd64eac..df6b77ea522 100644 --- a/include/ginkgo/core/solver/fcg.hpp +++ b/include/ginkgo/core/solver/fcg.hpp @@ -105,17 +105,10 @@ class Fcg void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Fcg(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Fcg(std::shared_ptr exec); explicit Fcg(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/gcr.hpp b/include/ginkgo/core/solver/gcr.hpp index 009a589a772..db34545fb00 100644 --- a/include/ginkgo/core/solver/gcr.hpp +++ b/include/ginkgo/core/solver/gcr.hpp @@ -117,21 +117,10 @@ class Gcr void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Gcr(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Gcr(std::shared_ptr exec); explicit Gcr(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - { - if (!parameters_.krylov_dim) { - parameters_.krylov_dim = gcr_default_krylov_dim; - } - } + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/gmres.hpp b/include/ginkgo/core/solver/gmres.hpp index 18fa37a9eff..acb782bf7f3 100644 --- a/include/ginkgo/core/solver/gmres.hpp +++ b/include/ginkgo/core/solver/gmres.hpp @@ -149,22 +149,10 @@ class Gmres void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit Gmres(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit Gmres(std::shared_ptr exec); explicit Gmres(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - { - if (!parameters_.krylov_dim) { - parameters_.krylov_dim = gmres_default_krylov_dim; - } - } + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/idr.hpp b/include/ginkgo/core/solver/idr.hpp index fc061e776a5..03efe95b81a 100644 --- a/include/ginkgo/core/solver/idr.hpp +++ b/include/ginkgo/core/solver/idr.hpp @@ -220,17 +220,10 @@ class Idr template void iterate(const VectorType* dense_b, VectorType* dense_x) const; - explicit Idr(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Idr(std::shared_ptr exec); explicit Idr(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/ir.hpp b/include/ginkgo/core/solver/ir.hpp index a76a3d247cb..1c3f8f33587 100644 --- a/include/ginkgo/core/solver/ir.hpp +++ b/include/ginkgo/core/solver/ir.hpp @@ -221,31 +221,10 @@ class Ir : public LinOp, void set_relaxation_factor( std::shared_ptr> new_factor); - explicit Ir(std::shared_ptr exec) : LinOp(std::move(exec)) - {} + explicit Ir(std::shared_ptr exec); explicit Ir(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnableSolverBase{std::move(system_matrix)}, - EnableIterativeBase{ - stop::combine(factory->get_parameters().criteria)}, - parameters_{factory->get_parameters()} - { - if (parameters_.generated_solver) { - this->set_solver(parameters_.generated_solver); - } else if (parameters_.solver) { - this->set_solver( - parameters_.solver->generate(this->get_system_matrix())); - } else { - this->set_solver(matrix::Identity::create( - this->get_executor(), this->get_size()[0])); - } - this->set_default_initial_guess(parameters_.default_initial_guess); - relaxation_factor_ = gko::initialize>( - {parameters_.relaxation_factor}, this->get_executor()); - } + std::shared_ptr system_matrix); private: std::shared_ptr solver_{}; diff --git a/include/ginkgo/core/solver/pipe_cg.hpp b/include/ginkgo/core/solver/pipe_cg.hpp index c271486766f..9dc1753bf61 100644 --- a/include/ginkgo/core/solver/pipe_cg.hpp +++ b/include/ginkgo/core/solver/pipe_cg.hpp @@ -114,18 +114,10 @@ class PipeCg void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override; - explicit PipeCg(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit PipeCg(std::shared_ptr exec); explicit PipeCg(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnablePreconditionedIterativeSolver>{ - std::move(system_matrix), factory->get_parameters()}, - parameters_{factory->get_parameters()} - {} + std::shared_ptr system_matrix); }; diff --git a/include/ginkgo/core/solver/triangular.hpp b/include/ginkgo/core/solver/triangular.hpp index 7256bb9c7ea..78a26d2e7ae 100644 --- a/include/ginkgo/core/solver/triangular.hpp +++ b/include/ginkgo/core/solver/triangular.hpp @@ -167,21 +167,10 @@ class LowerTrs : public LinOp, */ void generate(); - explicit LowerTrs(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit LowerTrs(std::shared_ptr exec); explicit LowerTrs(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnableSolverBase, CsrMatrix>{ - copy_and_convert_to(factory->get_executor(), - system_matrix)}, - parameters_{factory->get_parameters()} - { - this->generate(); - } + std::shared_ptr system_matrix); private: std::shared_ptr solve_struct_; @@ -334,21 +323,10 @@ class UpperTrs : public LinOp, */ void generate(); - explicit UpperTrs(std::shared_ptr exec) - : LinOp(std::move(exec)) - {} + explicit UpperTrs(std::shared_ptr exec); explicit UpperTrs(const Factory* factory, - std::shared_ptr system_matrix) - : LinOp(factory->get_executor(), - gko::transpose(system_matrix->get_size())), - EnableSolverBase, CsrMatrix>{ - copy_and_convert_to(factory->get_executor(), - system_matrix)}, - parameters_{factory->get_parameters()} - { - this->generate(); - } + std::shared_ptr system_matrix); private: std::shared_ptr solve_struct_;