From 384c764f6c4e13ac1c5d2bfda2c0828781036ca9 Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Wed, 19 Mar 2025 14:17:06 +0100 Subject: [PATCH 1/6] Allow to use LinOp as LSolverType in preconditioner::Ic if LSolverType is LinOp, value_type will be void and Ic will not give the default factory for l_solver and factorization --- core/config/preconditioner_ic_config.cpp | 8 +- core/preconditioner/ic.cpp | 21 ++- include/ginkgo/core/base/composition.hpp | 22 ++- include/ginkgo/core/preconditioner/ic.hpp | 191 +++++++++++++++------- reference/test/preconditioner/ic.cpp | 144 +++++++++++++++- 5 files changed, 311 insertions(+), 75 deletions(-) diff --git a/core/config/preconditioner_ic_config.cpp b/core/config/preconditioner_ic_config.cpp index 7d4cd3d9ca4..026faa5ec60 100644 --- a/core/config/preconditioner_ic_config.cpp +++ b/core/config/preconditioner_ic_config.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -97,6 +97,12 @@ deferred_factory_parameter parse( config, context, updated, make_type_selector(updated.get_value_typestr(), value_type_list()), make_type_selector(updated.get_index_typestr(), index_type_list())); + } else if (str == "LinOp") { + return dispatch::Configurator>( + config, context, updated, + // no effect + make_type_selector(updated.get_value_typestr(), value_type_list()), + make_type_selector(updated.get_index_typestr(), index_type_list())); } else { GKO_INVALID_CONFIG_VALUE("l_solver_type", str); } diff --git a/core/preconditioner/ic.cpp b/core/preconditioner/ic.cpp index 691795ad60b..665ba3e8fb6 100644 --- a/core/preconditioner/ic.cpp +++ b/core/preconditioner/ic.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -30,9 +30,15 @@ typename Ic::parameters_type ic_parse( auto params = Ic::build(); if (auto& obj = config.get("l_solver")) { - params.with_l_solver( - gko::config::parse_or_get_specific_factory< - const typename Ic::l_solver_type>(obj, context, td_for_child)); + if constexpr (std::is_same_v) { + params.with_l_solver( + gko::config::parse_or_get_factory( + obj, context, td_for_child)); + } else { + params.with_l_solver(gko::config::parse_or_get_specific_factory< + const typename Ic::l_solver_type>( + obj, context, td_for_child)); + } } if (auto& obj = config.get("factorization")) { params.with_factorization( @@ -73,6 +79,13 @@ GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GMRES_IC_PARSE); const config::type_descriptor&) GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_LOWERISAI_IC_PARSE); +#define GKO_DECLARE_LINOP_IC_PARSE(IndexType) \ + typename Ic::parameters_type \ + ic_parse>(const config::pnode&, \ + const config::registry&, \ + const config::type_descriptor&) +GKO_INSTANTIATE_FOR_EACH_INDEX_TYPE(GKO_DECLARE_LINOP_IC_PARSE); + } // namespace detail } // namespace preconditioner } // namespace gko diff --git a/include/ginkgo/core/base/composition.hpp b/include/ginkgo/core/base/composition.hpp index 9c16f8720aa..f10cf283986 100644 --- a/include/ginkgo/core/base/composition.hpp +++ b/include/ginkgo/core/base/composition.hpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -16,6 +16,21 @@ namespace gko { +/** + * CompositeionBase provides the interface + */ +class CompositionBase { +public: + /** + * Returns a list of operators of the composition. + * + * @return a list of operators + */ + virtual const std::vector>& get_operators() + const noexcept = 0; +}; + + /** * The Composition class can be used to compose linear operators `op1, op2, ..., * opn` and obtain the operator `op1 * op2 * ... * opn`. @@ -38,7 +53,8 @@ namespace gko { template class Composition : public EnableLinOp>, public EnableCreateMethod>, - public Transposable { + public Transposable, + public CompositionBase { friend class EnablePolymorphicObject; friend class EnableCreateMethod; @@ -52,7 +68,7 @@ class Composition : public EnableLinOp>, * @return a list of operators */ const std::vector>& get_operators() - const noexcept + const noexcept override { return operators_; } diff --git a/include/ginkgo/core/preconditioner/ic.hpp b/include/ginkgo/core/preconditioner/ic.hpp index 9260bfbb891..407e98c2d1e 100644 --- a/include/ginkgo/core/preconditioner/ic.hpp +++ b/include/ginkgo/core/preconditioner/ic.hpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -40,7 +40,8 @@ constexpr bool support_ic_parse = is_instantiation_of::value || is_instantiation_of::value || is_instantiation_of::value || - is_instantiation_of::value; + is_instantiation_of::value || + std::is_same_v; template < @@ -62,6 +63,54 @@ typename Ic::parameters_type ic_parse( const config::type_descriptor& td_for_child); +// helper for handle the transposed type of concrete type and LinOp +template +struct transposed_type_impl { + using type = typename Type::transposed_type; +}; + +template <> +struct transposed_type_impl { + using type = LinOp; +}; + + +template +using transposed_type = typename transposed_type_impl::type; + + +// helper to get factory type of concrete type or LinOp +template +struct factory_type_impl { + using type = typename Type::Factory; +}; + +template <> +struct factory_type_impl { + using type = LinOpFactory; +}; + + +template +using factory_type = typename factory_type_impl::type; + + +// helper to get value_type of concrete type or void for LinOp +template +struct get_value_type_impl { + using type = typename Type::value_type; +}; + +template <> +struct get_value_type_impl { + using type = void; +}; + + +template +using get_value_type = typename get_value_type_impl::type; + + } // namespace detail /** @@ -116,12 +165,13 @@ class Ic : public EnableLinOp>, public Transposable { public: static_assert( - std::is_same::value, + std::is_same< + detail::transposed_type>, + LSolverType>::value, "LSolverType::transposed_type must be symmetric"); - using value_type = typename LSolverType::value_type; + using value_type = detail::get_value_type; using l_solver_type = LSolverType; - using lh_solver_type = typename LSolverType::transposed_type; + using lh_solver_type = detail::transposed_type; using index_type = IndexType; using transposed_type = Ic; @@ -132,7 +182,7 @@ class Ic : public EnableLinOp>, public Transposable { /** * Factory for the L solver */ - std::shared_ptr + std::shared_ptr> l_solver_factory{}; /** @@ -142,14 +192,16 @@ class Ic : public EnableLinOp>, public Transposable { GKO_DEPRECATED("use with_l_solver instead") parameters_type& with_l_solver_factory( - deferred_factory_parameter + deferred_factory_parameter< + const detail::factory_type> solver) { return with_l_solver(std::move(solver)); } parameters_type& with_l_solver( - deferred_factory_parameter + deferred_factory_parameter< + const detail::factory_type> solver) { this->l_solver_generator = std::move(solver); @@ -185,7 +237,7 @@ class Ic : public EnableLinOp>, public Transposable { } private: - deferred_factory_parameter + deferred_factory_parameter> l_solver_generator; deferred_factory_parameter factorization_generator; @@ -204,16 +256,21 @@ class Ic : public EnableLinOp>, public Transposable { * @param context the registry * @param td_for_child the type descriptor for children configs. The * default uses the value/index type of this class. + * When l_solver_type uses LinOp not concrete type, it + * will use the default_precision in ginkgo. * * @return parameters * * @note only support the following for l_solver: - * Ir, Gmres, LowerTrs, and LowerIsai + * Ir, Gmres, LowerTrs, LowerIsai, LinOp */ static parameters_type parse( const config::pnode& config, const config::registry& context, const config::type_descriptor& td_for_child = - config::make_type_descriptor()) + config::make_type_descriptor< + std::conditional_t, + default_precision, value_type>, + index_type>()) { return detail::ic_parse(config, context, td_for_child); } @@ -244,11 +301,11 @@ class Ic : public EnableLinOp>, public Transposable { new transposed_type{this->get_executor()}}; transposed->set_size(gko::transpose(this->get_size())); transposed->l_solver_ = - share(as( - this->get_lh_solver()->transpose())); + share(as>( + as(this->get_lh_solver())->transpose())); transposed->lh_solver_ = - share(as( - this->get_l_solver()->transpose())); + share(as>( + as(this->get_l_solver())->transpose())); return std::move(transposed); } @@ -259,11 +316,11 @@ class Ic : public EnableLinOp>, public Transposable { new transposed_type{this->get_executor()}}; transposed->set_size(gko::transpose(this->get_size())); transposed->l_solver_ = - share(as( - this->get_lh_solver()->conj_transpose())); + share(as>( + as(this->get_lh_solver())->conj_transpose())); transposed->lh_solver_ = - share(as( - this->get_l_solver()->conj_transpose())); + share(as>( + as(this->get_l_solver())->conj_transpose())); return std::move(transposed); } @@ -327,30 +384,22 @@ class Ic : public EnableLinOp>, public Transposable { protected: void apply_impl(const LinOp* b, LinOp* x) const override { - // take care of real-to-complex apply - precision_dispatch_real_complex( - [&](auto dense_b, auto dense_x) { - this->set_cache_to(dense_b); - l_solver_->apply(dense_b, cache_.intermediate); - if (lh_solver_->apply_uses_initial_guess()) { - dense_x->copy_from(cache_.intermediate); - } - lh_solver_->apply(cache_.intermediate, dense_x); - }, - b, x); + // let the solver to handle the precision + this->set_cache_to(b); + l_solver_->apply(b, cache_.intermediate); + if (lh_solver_->apply_uses_initial_guess()) { + x->copy_from(cache_.intermediate); + } + lh_solver_->apply(cache_.intermediate, x); } void apply_impl(const LinOp* alpha, const LinOp* b, const LinOp* beta, LinOp* x) const override { - precision_dispatch_real_complex( - [&](auto dense_alpha, auto dense_b, auto dense_beta, auto dense_x) { - this->set_cache_to(dense_b); - l_solver_->apply(dense_b, cache_.intermediate); - lh_solver_->apply(dense_alpha, cache_.intermediate, dense_beta, - dense_x); - }, - alpha, b, beta, x); + // let the solver to handle the precision + this->set_cache_to(b); + l_solver_->apply(b, cache_.intermediate); + lh_solver_->apply(alpha, cache_.intermediate, beta, x); } explicit Ic(std::shared_ptr exec) @@ -361,24 +410,28 @@ class Ic : public EnableLinOp>, public Transposable { : EnableLinOp(factory->get_executor(), lin_op->get_size()), parameters_{factory->get_parameters()} { - auto comp = - std::dynamic_pointer_cast>(lin_op); + auto comp = std::dynamic_pointer_cast(lin_op); std::shared_ptr l_factor; // build factorization if we weren't passed a composition if (!comp) { auto exec = lin_op->get_executor(); if (!parameters_.factorization_factory) { - parameters_.factorization_factory = - factorization::ParIc::build() - .with_both_factors(false) - .on(exec); + if constexpr (std::is_same_v) { + // linop does not have the value_type such that we do not + // know which par_ic can use. + GKO_NOT_IMPLEMENTED; + } else { + parameters_.factorization_factory = + factorization::ParIc::build() + .with_both_factors(false) + .on(exec); + } } auto fact = std::shared_ptr( parameters_.factorization_factory->generate(lin_op)); // ensure that the result is a composition - comp = - std::dynamic_pointer_cast>(fact); + comp = std::dynamic_pointer_cast(fact); if (!comp) { GKO_NOT_SUPPORTED(comp); } @@ -394,19 +447,28 @@ class Ic : public EnableLinOp>, public Transposable { // If no factories are provided, generate default ones if (!parameters_.l_solver_factory) { - l_solver_ = generate_default_solver(exec, l_factor); - // If comp contains both factors: use the transposed factor to avoid - // transposing twice - if (comp->get_operators().size() == 2) { - auto lh_factor = comp->get_operators()[1]; - GKO_ASSERT_EQUAL_DIMENSIONS(l_factor, lh_factor); - lh_solver_ = as(l_solver_->conj_transpose()); + if constexpr (std::is_same_v) { + // linop can not have a default factory generation + GKO_NOT_IMPLEMENTED; } else { - lh_solver_ = as(l_solver_->conj_transpose()); + l_solver_ = + generate_default_solver(exec, l_factor); + // If comp contains both factors: use the transposed factor to + // avoid transposing twice + if (comp->get_operators().size() == 2) { + auto lh_factor = comp->get_operators()[1]; + GKO_ASSERT_EQUAL_DIMENSIONS(l_factor, lh_factor); + lh_solver_ = as( + as(l_solver_)->conj_transpose()); + } else { + lh_solver_ = as( + as(l_solver_)->conj_transpose()); + } } } else { l_solver_ = parameters_.l_solver_factory->generate(l_factor); - lh_solver_ = as(l_solver_->conj_transpose()); + lh_solver_ = as( + as(l_solver_)->conj_transpose()); } } @@ -419,12 +481,13 @@ class Ic : public EnableLinOp>, public Transposable { */ void set_cache_to(const LinOp* b) const { - if (cache_.intermediate == nullptr) { - cache_.intermediate = - matrix::Dense::create(this->get_executor()); + if (cache_.intermediate == nullptr || + cache_.intermediate->get_size() != b->get_size()) { + cache_.intermediate = b->clone(this->get_executor()); + } else { + // Use b as the initial guess for the first triangular solve + cache_.intermediate->copy_from(b); } - // Use b as the initial guess for the first triangular solve - cache_.intermediate->copy_from(b); } @@ -436,7 +499,8 @@ class Ic : public EnableLinOp>, public Transposable { * preconditioner. */ template - static std::enable_if_t::value, + static std::enable_if_t::value && + !std::is_same_v, std::unique_ptr> generate_default_solver(const std::shared_ptr& exec, const std::shared_ptr& mtx) @@ -458,7 +522,8 @@ class Ic : public EnableLinOp>, public Transposable { * @copydoc generate_default_solver */ template - static std::enable_if_t::value, + static std::enable_if_t::value && + !std::is_same_v, std::unique_ptr> generate_default_solver(const std::shared_ptr& exec, const std::shared_ptr& mtx) diff --git a/reference/test/preconditioner/ic.cpp b/reference/test/preconditioner/ic.cpp index 16ffc8d7b3c..3858c2dda5f 100644 --- a/reference/test/preconditioner/ic.cpp +++ b/reference/test/preconditioner/ic.cpp @@ -1,4 +1,4 @@ -// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// SPDX-FileCopyrightText: 2017 - 2025 The Ginkgo authors // // SPDX-License-Identifier: BSD-3-Clause @@ -30,11 +30,11 @@ class Ic : public ::testing::Test { typename std::tuple_element<1, decltype(ValueIndexType())>::type; using Mtx = gko::matrix::Csr; using Vec = gko::matrix::Dense; - using ic_prec_type = - gko::preconditioner::Ic, - index_type>; + using lower_trs = gko::solver::LowerTrs; + using ic_prec_type = gko::preconditioner::Ic; using ic_isai_prec_type = gko::preconditioner::Ic< gko::preconditioner::LowerIsai, index_type>; + using ic_linop_prec_type = gko::preconditioner::Ic; using Composition = gko::Composition; Ic() @@ -51,6 +51,13 @@ class Ic : public ::testing::Test { l_lh_composition(Composition::create(l_factor, lh_factor)), ic_pre_factory(ic_prec_type::build().on(exec)), ic_isai_pre_factory(ic_isai_prec_type::build().on(exec)), + ic_linop_pre_factory( + ic_linop_prec_type::build() + .with_l_solver(lower_trs::build()) + .with_factorization( + gko::factorization::ParIc::build() + .with_both_factors(false)) + .on(exec)), tol{r::value} {} @@ -64,6 +71,7 @@ class Ic : public ::testing::Test { std::shared_ptr l_lh_composition; std::shared_ptr ic_pre_factory; std::shared_ptr ic_isai_pre_factory; + std::shared_ptr ic_linop_pre_factory; gko::remove_complex tol; }; @@ -126,6 +134,72 @@ TYPED_TEST(Ic, BuildsIsaiFromMatrix) } +TYPED_TEST(Ic, BuildsTwoFactorCompositionWithLinOp) +{ + using lower_trs = typename TestFixture::lower_trs; + + auto precond = this->ic_linop_pre_factory->generate(this->l_lh_composition); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver()) + ->get_system_matrix(), + this->lh_factor, this->tol); +} + + +TYPED_TEST(Ic, BuildsOneFactorCompositionWithLinOp) +{ + using lower_trs = typename TestFixture::lower_trs; + + auto precond = this->ic_linop_pre_factory->generate(this->l_composition); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver()) + ->get_system_matrix(), + this->lh_factor, this->tol); +} + + +TYPED_TEST(Ic, BuildsSymmetricTwoFactorWithLinOp) +{ + using lower_trs = typename TestFixture::lower_trs; + + auto precond = this->ic_linop_pre_factory->generate(this->l_lh_composition); + + // the first factor should be identical, the second not as it was transposed + ASSERT_EQ( + gko::as(precond->get_l_solver())->get_system_matrix().get(), + this->l_factor.get()); + ASSERT_NE( + gko::as(precond->get_lh_solver()) + ->get_system_matrix() + .get(), + this->lh_factor.get()); +} + + +TYPED_TEST(Ic, BuildsFromMatrixWithLinOp) +{ + using lower_trs = typename TestFixture::lower_trs; + + auto precond = this->ic_linop_pre_factory->generate(this->mtx); + + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_l_solver())->get_system_matrix(), + this->l_factor, this->tol); + GKO_ASSERT_MTX_NEAR( + gko::as(precond->get_lh_solver()) + ->get_system_matrix(), + this->lh_factor, this->tol); +} + + TYPED_TEST(Ic, ThrowOnWrongCompositionInput) { using Composition = typename TestFixture::Composition; @@ -226,6 +300,52 @@ TYPED_TEST(Ic, CanBeConjTransposed) } +TYPED_TEST(Ic, CanBeTransposedWithLinOp) +{ + using Ic = typename TestFixture::ic_linop_prec_type; + using Mtx = typename TestFixture::Mtx; + using lower_trs = typename TestFixture::lower_trs; + using upper_trs = typename lower_trs::transposed_type; + auto ic = this->ic_linop_pre_factory->generate(this->l_lh_composition); + auto l_ref = gko::as( + gko::as(ic->get_l_solver())->get_system_matrix()); + auto lh_ref = gko::as( + gko::as(ic->get_lh_solver())->get_system_matrix()); + + auto transp = gko::as(ic->transpose()); + + auto l_transp = gko::as( + gko::as(transp->get_l_solver())->get_system_matrix()); + auto lh_transp = gko::as( + gko::as(transp->get_lh_solver())->get_system_matrix()); + GKO_ASSERT_MTX_NEAR(l_transp, l_ref, 0); + GKO_ASSERT_MTX_NEAR(lh_transp, lh_ref, 0); +} + + +TYPED_TEST(Ic, CanBeConjTransposedWithLinOp) +{ + using Ic = typename TestFixture::ic_linop_prec_type; + using Mtx = typename TestFixture::Mtx; + using lower_trs = typename TestFixture::lower_trs; + using upper_trs = typename lower_trs::transposed_type; + auto ic = this->ic_linop_pre_factory->generate(this->l_lh_composition); + auto l_ref = gko::as( + gko::as(ic->get_l_solver())->get_system_matrix()); + auto lh_ref = gko::as( + gko::as(ic->get_lh_solver())->get_system_matrix()); + + auto transp = gko::as(ic->conj_transpose()); + + auto l_transp = gko::as( + gko::as(transp->get_l_solver())->get_system_matrix()); + auto lh_transp = gko::as( + gko::as(transp->get_lh_solver())->get_system_matrix()); + GKO_ASSERT_MTX_NEAR(l_transp, l_ref, 0); + GKO_ASSERT_MTX_NEAR(lh_transp, lh_ref, 0); +} + + TYPED_TEST(Ic, SolvesSingleRhs) { using ic_prec_type = typename TestFixture::ic_prec_type; @@ -392,4 +512,20 @@ TYPED_TEST(Ic, SolvesMultipleRhs) } +TYPED_TEST(Ic, SolvesMultipleRhsWithLinop) +{ + using Vec = typename TestFixture::Vec; + const auto b = gko::initialize( + {{1.0, 2.0, 3.0}, {3.0, 6.0, 9.0}, {6.0, 12.0, 18.0}}, this->exec); + auto x = Vec::create(this->exec, gko::dim<2>{3, 3}); + auto preconditioner = this->ic_linop_pre_factory->generate(this->mtx); + + preconditioner->apply(b, x); + + GKO_ASSERT_MTX_NEAR( + x, l({{3.0, 6.0, 9.0}, {-2.0, -4.0, -6.0}, {4.0, 8.0, 12.0}}), + this->tol); +} + + } // namespace From 70e584cec495be003ce1a165b9787f5cdb99368f Mon Sep 17 00:00:00 2001 From: "Yu-Hsiang M. Tsai" Date: Wed, 19 Mar 2025 16:10:02 +0100 Subject: [PATCH 2/6] only allow LinOp in IC during parse --- core/config/preconditioner_ic_config.cpp | 77 +++-------------------- core/preconditioner/ic.cpp | 30 --------- core/test/config/preconditioner.cpp | 15 +++-- include/ginkgo/core/preconditioner/ic.hpp | 13 ++-- 4 files changed, 21 insertions(+), 114 deletions(-) diff --git a/core/config/preconditioner_ic_config.cpp b/core/config/preconditioner_ic_config.cpp index 026faa5ec60..0900bcb7d6e 100644 --- a/core/config/preconditioner_ic_config.cpp +++ b/core/config/preconditioner_ic_config.cpp @@ -22,90 +22,31 @@ namespace gko { namespace config { -// For Ic and Ilu, we use additional ValueType to help Solver type decision -template class IcSolverHelper { public: - template + template class Configurator { public: - static - typename gko::preconditioner::Ic::parameters_type - parse(const pnode& config, const registry& context, - const type_descriptor& td_for_child) + static typename gko::preconditioner::Ic::parameters_type + parse(const pnode& config, const registry& context, + const type_descriptor& td_for_child) { - return gko::preconditioner::Ic::parse( + return gko::preconditioner::Ic::parse( config, context, td_for_child); } }; }; -// Do not use the partial specialization for SolverBase and SolverBase -// because the default template arguments are allowed for a template template -// argument (detail: CWG 150 after c++17 -// https://en.cppreference.com/w/cpp/language/template_parameters#Template_template_arguments) -template