From de40cea082e67e06f46b699cfb79874131c30ba3 Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Thu, 11 Jul 2024 07:24:09 +0000 Subject: [PATCH] [sor] review updates: - documentation - tests - don't build upper solver if not symmetric - remove unused includes - use identity function Co-authored-by: Yu-Hsiang M. Tsai Co-authored-by: Fritz Goebel --- core/preconditioner/sor.cpp | 5 +- core/test/config/preconditioner.cpp | 65 ++++++++++++++----- .../core/preconditioner/gauss_seidel.hpp | 18 +++-- include/ginkgo/core/preconditioner/sor.hpp | 19 ++++-- reference/CMakeLists.txt | 2 +- reference/preconditioner/sor_kernels.cpp | 2 +- reference/test/preconditioner/sor_kernels.cpp | 32 ++++----- 7 files changed, 98 insertions(+), 45 deletions(-) diff --git a/core/preconditioner/sor.cpp b/core/preconditioner/sor.cpp index 0ff534a268c..c9905c5f73c 100644 --- a/core/preconditioner/sor.cpp +++ b/core/preconditioner/sor.cpp @@ -96,10 +96,11 @@ std::unique_ptr Sor::generate_impl( auto l_trs_factory = parameters_.l_solver ? parameters_.l_solver : LTrs::build().on(exec); - auto u_trs_factory = - parameters_.u_solver ? parameters_.u_solver : UTrs::build().on(exec); if (parameters_.symmetric) { + auto u_trs_factory = parameters_.u_solver ? parameters_.u_solver + : UTrs::build().on(exec); + array l_row_ptrs{exec, size[0] + 1}; array u_row_ptrs{exec, size[0] + 1}; exec->run(make_initialize_row_ptrs_l_u( diff --git a/core/test/config/preconditioner.cpp b/core/test/config/preconditioner.cpp index c3941f504e2..bbb6ff2ffd7 100644 --- a/core/test/config/preconditioner.cpp +++ b/core/test/config/preconditioner.cpp @@ -36,6 +36,7 @@ struct PreconditionerConfigTest { using changed_type = ChangedType; using default_type = DefaultType; using preconditioner_config_test = PreconditionerConfigTest; + static pnode::map_type setup_base() { return {{"type", pnode{"preconditioner::Ic"}}}; @@ -328,12 +329,23 @@ struct Sor config_map["relaxation_factor"] = pnode{0.8}; // float can be cast to double without issues param.with_relaxation_factor(0.8f); - config_map["l_solver"] = pnode{ - {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_l_solver(DummyIr::build()); - config_map["u_solver"] = pnode{ - {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_u_solver(DummyIr::build()); + if (from_reg) { + config_map["l_solver"] = pnode{"l_solver"}; + param.with_l_solver( + detail::registry_accessor::get_data( + reg, "l_solver")); + config_map["u_solver"] = pnode{"u_solver"}; + param.with_u_solver( + detail::registry_accessor::get_data( + reg, "u_solver")); + } else { + config_map["l_solver"] = pnode{{{"type", pnode{"solver::Ir"}}, + {"value_type", pnode{"float32"}}}}; + param.with_l_solver(Ir::build()); + config_map["u_solver"] = pnode{{{"type", pnode{"solver::Ir"}}, + {"value_type", pnode{"float32"}}}}; + param.with_u_solver(Ir::build()); + } } template @@ -345,8 +357,13 @@ struct Sor ASSERT_EQ(res_param.skip_sorting, ans_param.skip_sorting); ASSERT_EQ(res_param.symmetric, ans_param.symmetric); ASSERT_EQ(res_param.relaxation_factor, ans_param.relaxation_factor); - ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver)); - ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver)); + if (from_reg) { + ASSERT_EQ(res_param.l_solver, ans_param.l_solver); + ASSERT_EQ(res_param.u_solver, ans_param.u_solver); + } else { + ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver)); + ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver)); + } } }; @@ -375,12 +392,23 @@ struct GaussSeidel param.with_skip_sorting(true); config_map["symmetric"] = pnode{true}; param.with_symmetric(true); - config_map["l_solver"] = pnode{ - {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_l_solver(DummyIr::build()); - config_map["u_solver"] = pnode{ - {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_u_solver(DummyIr::build()); + if (from_reg) { + config_map["l_solver"] = pnode{"l_solver"}; + param.with_l_solver( + detail::registry_accessor::get_data( + reg, "l_solver")); + config_map["u_solver"] = pnode{"u_solver"}; + param.with_u_solver( + detail::registry_accessor::get_data( + reg, "u_solver")); + } else { + config_map["l_solver"] = pnode{{{"type", pnode{"solver::Ir"}}, + {"value_type", pnode{"float32"}}}}; + param.with_l_solver(Ir::build()); + config_map["u_solver"] = pnode{{{"type", pnode{"solver::Ir"}}, + {"value_type", pnode{"float32"}}}}; + param.with_u_solver(Ir::build()); + } } template @@ -391,8 +419,13 @@ struct GaussSeidel ASSERT_EQ(res_param.skip_sorting, ans_param.skip_sorting); ASSERT_EQ(res_param.symmetric, ans_param.symmetric); - ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver)); - ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver)); + if (from_reg) { + ASSERT_EQ(res_param.l_solver, ans_param.l_solver); + ASSERT_EQ(res_param.u_solver, ans_param.u_solver); + } else { + ASSERT_EQ(typeid(res_param.l_solver), typeid(ans_param.l_solver)); + ASSERT_EQ(typeid(res_param.u_solver), typeid(ans_param.u_solver)); + } } }; diff --git a/include/ginkgo/core/preconditioner/gauss_seidel.hpp b/include/ginkgo/core/preconditioner/gauss_seidel.hpp index 7003dd54740..38ae0a260ea 100644 --- a/include/ginkgo/core/preconditioner/gauss_seidel.hpp +++ b/include/ginkgo/core/preconditioner/gauss_seidel.hpp @@ -6,8 +6,6 @@ #define GKO_PUBLIC_CORE_PRECONDITIONER_GAUSS_SEIDEL_HPP_ -#include - #include #include #include @@ -22,7 +20,8 @@ namespace preconditioner { /** * This class generates the Gauss-Seidel preconditioner. * - * This is the special case of $\omega = 1$ of the (S)SOR preconditioner. + * This is the special case of the relaxation factor $\omega = 1$ of the (S)SOR + * preconditioner. * * @see Sor */ @@ -49,12 +48,12 @@ class GaussSeidel // determines if Gauss-Seidel or symmetric Gauss-Seidel should be used bool GKO_FACTORY_PARAMETER_SCALAR(symmetric, false); - // factory for the lower triangular factor solver + // factory for the lower triangular factor solver, defaults to LowerTrs std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( l_solver); // factory for the upper triangular factor solver, unused if symmetric - // is false + // is false, defaults to UpperTrs std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( u_solver); }; @@ -84,6 +83,15 @@ class GaussSeidel /** Creates a new parameter_type to set up the factory. */ static parameters_type build() { return {}; } + /** + * Creates a new parameter_type based on an input configuration. + * + * @param config The parameter tree for the Gauss-Seidel. + * @param context The context for the parsing. + * @param td_for_child The type_descriptor used for child configurations. + * + * @return a parameter_type object with the settings from config + */ static parameters_type parse( const config::pnode& config, const config::registry& context, const config::type_descriptor& td_for_child = diff --git a/include/ginkgo/core/preconditioner/sor.hpp b/include/ginkgo/core/preconditioner/sor.hpp index 276d718dacb..c35dc631e61 100644 --- a/include/ginkgo/core/preconditioner/sor.hpp +++ b/include/ginkgo/core/preconditioner/sor.hpp @@ -6,8 +6,6 @@ #define GKO_PUBLIC_CORE_PRECONDITIONER_SOR_HPP_ -#include - #include #include #include @@ -36,6 +34,8 @@ namespace preconditioner { * M = \frac{1}{\omega (2 - \omega)} (D + \omega L) D^{-1} (D + \omega U) , * \quad 0 < \omega < 2. * $$ + * A detailed description can be found in Iterative Methods for Sparse Linear + * Systems (Y. Saad) ch. 4.1. * * This class is a factory, which will only generate the preconditioner. The * resulting LinOp will represent the application of $M^{-1}$. @@ -69,12 +69,12 @@ class Sor remove_complex GKO_FACTORY_PARAMETER_SCALAR( relaxation_factor, remove_complex(1.2)); - // factory for the lower triangular factor solver + // factory for the lower triangular factor solver, defaults to LowerTrs std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( l_solver); // factory for the upper triangular factor solver, unused if symmetric - // is false + // is false, defaults to UpperTrs std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( u_solver); }; @@ -104,6 +104,15 @@ class Sor /** Creates a new parameter_type to set up the factory. */ static parameters_type build() { return {}; } + /** + * Creates a new parameter_type based on an input configuration. + * + * @param config The parameter tree for the SOR. + * @param context The context for the parsing. + * @param td_for_child The type_descriptor used for child configurations. + * + * @return a parameter_type object with the settings from config + */ static parameters_type parse( const config::pnode& config, const config::registry& context, const config::type_descriptor& td_for_child = @@ -124,6 +133,8 @@ class Sor private: parameters_type parameters_; }; + + } // namespace preconditioner } // namespace gko diff --git a/reference/CMakeLists.txt b/reference/CMakeLists.txt index ab6c210518b..e2f27dab57e 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -43,7 +43,7 @@ target_sources(ginkgo_reference matrix/sparsity_csr_kernels.cpp multigrid/pgm_kernels.cpp preconditioner/batch_jacobi_kernels.cpp - preconditioner/sor_kernels.cpp + preconditioner/sor_kernels.cpp preconditioner/isai_kernels.cpp preconditioner/jacobi_kernels.cpp reorder/rcm_kernels.cpp diff --git a/reference/preconditioner/sor_kernels.cpp b/reference/preconditioner/sor_kernels.cpp index 88ac422dd02..9c109788318 100644 --- a/reference/preconditioner/sor_kernels.cpp +++ b/reference/preconditioner/sor_kernels.cpp @@ -29,7 +29,7 @@ void initialize_weighted_l( system_matrix, l_mtx, factorization::helpers::triangular_mtx_closure( [inv_weight](auto val) { return val * inv_weight; }, - [](auto val) { return val; })); + factorization::helpers::identity{})); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( diff --git a/reference/test/preconditioner/sor_kernels.cpp b/reference/test/preconditioner/sor_kernels.cpp index f2bf5f186f9..18c055aa6d9 100644 --- a/reference/test/preconditioner/sor_kernels.cpp +++ b/reference/test/preconditioner/sor_kernels.cpp @@ -80,15 +80,15 @@ TYPED_TEST(Sor, CanInitializeLFactorWithWeight) result->scale( gko::initialize>({0.0}, this->exec)); std::shared_ptr expected_l = - gko::initialize({{1, 0, 0, 0, 0}, - {-2, 1, 0, 0, 0}, - {0, -5, 1, 0, 0}, - {-3, 0, 0, 1, 0}, - {-4, 0, -6, -7, 1}}, + gko::initialize({{2 * this->diag_value, 0, 0, 0, 0}, + {-2, 2 * this->diag_value, 0, 0, 0}, + {0, -5, 2 * this->diag_value, 0, 0}, + {-3, 0, 0, 2 * this->diag_value, 0}, + {-4, 0, -6, -7, 2 * this->diag_value}}, this->exec); gko::kernels::reference::sor::initialize_weighted_l( - this->exec, this->mtx.get(), this->diag_value, result.get()); + this->exec, this->mtx.get(), 0.5f, result.get()); GKO_ASSERT_MTX_NEAR(result, expected_l, r::value); } @@ -122,15 +122,16 @@ TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight) gko::initialize>({0.0}, this->exec)); result_u->scale( gko::initialize>({0.0}, this->exec)); - auto diag_weight = static_cast>( - 1.0 / (2 - this->diag_value)); - auto off_diag_weight = this->diag_value * diag_weight; + auto factor = static_cast>(0.5); + auto diag_weight = + static_cast>(1.0 / (2 - factor)); + auto off_diag_weight = factor * diag_weight; std::shared_ptr expected_l = - gko::initialize({{1, 0, 0, 0, 0}, - {-2, 1, 0, 0, 0}, - {0, -5, 1, 0, 0}, - {-3, 0, 0, 1, 0}, - {-4, 0, -6, -7, 1}}, + gko::initialize({{2 * this->diag_value, 0, 0, 0, 0}, + {-2, 2 * this->diag_value, 0, 0, 0}, + {0, -5, 2 * this->diag_value, 0, 0}, + {-3, 0, 0, 2 * this->diag_value, 0}, + {-4, 0, -6, -7, 2 * this->diag_value}}, this->exec); std::shared_ptr expected_u = gko::initialize( {{this->diag_value * diag_weight, 2 * off_diag_weight, 0, @@ -142,8 +143,7 @@ TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight) this->exec); gko::kernels::reference::sor::initialize_weighted_l_u( - this->exec, this->mtx.get(), this->diag_value, result_l.get(), - result_u.get()); + this->exec, this->mtx.get(), factor, result_l.get(), result_u.get()); GKO_ASSERT_MTX_NEAR(result_l, expected_l, r::value); GKO_ASSERT_MTX_NEAR(result_u, expected_u, r::value);