From 6d5a3f428964f90a82da9b889475d9e2b1f95307 Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Tue, 25 Jun 2024 13:46:55 +0200 Subject: [PATCH 1/7] [fact] extract l/u initialization --- core/factorization/factorization_helpers.hpp | 56 +++++++++ .../factorization/factorization_helpers.hpp | 112 ++++++++++++++++++ .../factorization/factorization_kernels.cpp | 97 +++------------ 3 files changed, 187 insertions(+), 78 deletions(-) create mode 100644 core/factorization/factorization_helpers.hpp create mode 100644 reference/factorization/factorization_helpers.hpp diff --git a/core/factorization/factorization_helpers.hpp b/core/factorization/factorization_helpers.hpp new file mode 100644 index 00000000000..16ead3a198d --- /dev/null +++ b/core/factorization/factorization_helpers.hpp @@ -0,0 +1,56 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GINKGO_CORE_FACTORIZATION_FACTORIZATION_HELPERS_HPP +#define GINKGO_CORE_FACTORIZATION_FACTORIZATION_HELPERS_HPP + + +#include + + +namespace gko { +namespace factorization { + + +struct identity { + template + constexpr T operator()(T value) + { + return value; + } +}; + + +template +class triangular_mtx_closure { +public: + constexpr triangular_mtx_closure(DiagClosure diag_closure, + OffDiagClosure off_diag_closure) + : diag_closure_(std::move(diag_closure)), + off_diag_closure_(std::move(off_diag_closure)) + {} + + template + constexpr T map_diag(T value) + { + return diag_closure_(value); + } + + template + constexpr T map_off_diag(T value) + { + return off_diag_closure_(value); + } + +private: + DiagClosure diag_closure_; + OffDiagClosure off_diag_closure_; +}; + + +} // namespace factorization +} // namespace gko + + +#endif // GINKGO_CORE_FACTORIZATION_FACTORIZATION_HELPERS_HPP diff --git a/reference/factorization/factorization_helpers.hpp b/reference/factorization/factorization_helpers.hpp new file mode 100644 index 00000000000..9b5afb518af --- /dev/null +++ b/reference/factorization/factorization_helpers.hpp @@ -0,0 +1,112 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include "core/factorization/factorization_helpers.hpp" + + +namespace gko { +namespace kernels { +namespace reference { +namespace factorization { +namespace helpers { + + +using namespace ::gko::factorization; + + +template +void initialize_l_u(const matrix::Csr* system_matrix, + matrix::Csr* csr_l, + matrix::Csr* csr_u, + LClosure&& l_closure, UClosure&& u_closure) +{ + const auto row_ptrs = system_matrix->get_const_row_ptrs(); + const auto col_idxs = system_matrix->get_const_col_idxs(); + const auto vals = system_matrix->get_const_values(); + + const auto row_ptrs_l = csr_l->get_const_row_ptrs(); + auto col_idxs_l = csr_l->get_col_idxs(); + auto vals_l = csr_l->get_values(); + + const auto row_ptrs_u = csr_u->get_const_row_ptrs(); + auto col_idxs_u = csr_u->get_col_idxs(); + auto vals_u = csr_u->get_values(); + + for (size_type row = 0; row < system_matrix->get_size()[0]; ++row) { + size_type current_index_l = row_ptrs_l[row]; + size_type current_index_u = + row_ptrs_u[row] + 1; // we treat the diagonal separately + // if there is no diagonal value, set it to 1 by default + auto diag_val = one(); + for (size_type el = row_ptrs[row]; el < row_ptrs[row + 1]; ++el) { + const auto col = col_idxs[el]; + const auto val = vals[el]; + if (col < row) { + col_idxs_l[current_index_l] = col; + vals_l[current_index_l] = l_closure.map_off_diag(val); + ++current_index_l; + } else if (col == row) { + // save diagonal value + diag_val = val; + } else { // col > row + col_idxs_u[current_index_u] = col; + vals_u[current_index_u] = u_closure.map_off_diag(val); + ++current_index_u; + } + } + // store diagonal values separately + auto l_diag_idx = row_ptrs_l[row + 1] - 1; + auto u_diag_idx = row_ptrs_u[row]; + col_idxs_l[l_diag_idx] = row; + col_idxs_u[u_diag_idx] = row; + vals_l[l_diag_idx] = l_closure.map_diag(diag_val); + vals_u[u_diag_idx] = u_closure.map_diag(diag_val); + } +} + + +template +void initialize_l(const matrix::Csr* system_matrix, + matrix::Csr* csr_l, Closure&& closure) +{ + const auto row_ptrs = system_matrix->get_const_row_ptrs(); + const auto col_idxs = system_matrix->get_const_col_idxs(); + const auto vals = system_matrix->get_const_values(); + + const auto row_ptrs_l = csr_l->get_const_row_ptrs(); + auto col_idxs_l = csr_l->get_col_idxs(); + auto vals_l = csr_l->get_values(); + + for (size_type row = 0; row < system_matrix->get_size()[0]; ++row) { + size_type current_index_l = row_ptrs_l[row]; + // if there is no diagonal value, set it to 1 by default + auto diag_val = one(); + for (size_type el = row_ptrs[row]; el < row_ptrs[row + 1]; ++el) { + const auto col = col_idxs[el]; + const auto val = vals[el]; + if (col < row) { + col_idxs_l[current_index_l] = col; + vals_l[current_index_l] = closure.map_off_diag(val); + ++current_index_l; + } else if (col == row) { + // save diagonal value + diag_val = val; + } + } + // store diagonal values separately + auto l_diag_idx = row_ptrs_l[row + 1] - 1; + col_idxs_l[l_diag_idx] = row; + vals_l[l_diag_idx] = closure.map_diag(diag_val); + } +} + + +} // namespace helpers +} // namespace factorization +} // namespace reference +} // namespace kernels +} // namespace gko diff --git a/reference/factorization/factorization_kernels.cpp b/reference/factorization/factorization_kernels.cpp index 085e2f62ecc..99b522ffba9 100644 --- a/reference/factorization/factorization_kernels.cpp +++ b/reference/factorization/factorization_kernels.cpp @@ -12,6 +12,7 @@ #include "core/components/prefix_sum_kernels.hpp" #include "core/matrix/csr_builder.hpp" +#include "reference/factorization/factorization_helpers.hpp" namespace gko { @@ -168,48 +169,12 @@ void initialize_l_u(std::shared_ptr exec, matrix::Csr* csr_l, matrix::Csr* csr_u) { - const auto row_ptrs = system_matrix->get_const_row_ptrs(); - const auto col_idxs = system_matrix->get_const_col_idxs(); - const auto vals = system_matrix->get_const_values(); - - const auto row_ptrs_l = csr_l->get_const_row_ptrs(); - auto col_idxs_l = csr_l->get_col_idxs(); - auto vals_l = csr_l->get_values(); - - const auto row_ptrs_u = csr_u->get_const_row_ptrs(); - auto col_idxs_u = csr_u->get_col_idxs(); - auto vals_u = csr_u->get_values(); - - for (size_type row = 0; row < system_matrix->get_size()[0]; ++row) { - size_type current_index_l = row_ptrs_l[row]; - size_type current_index_u = - row_ptrs_u[row] + 1; // we treat the diagonal separately - // if there is no diagonal value, set it to 1 by default - auto diag_val = one(); - for (size_type el = row_ptrs[row]; el < row_ptrs[row + 1]; ++el) { - const auto col = col_idxs[el]; - const auto val = vals[el]; - if (col < row) { - col_idxs_l[current_index_l] = col; - vals_l[current_index_l] = val; - ++current_index_l; - } else if (col == row) { - // save diagonal value - diag_val = val; - } else { // col > row - col_idxs_u[current_index_u] = col; - vals_u[current_index_u] = val; - ++current_index_u; - } - } - // store diagonal values separately - auto l_diag_idx = row_ptrs_l[row + 1] - 1; - auto u_diag_idx = row_ptrs_u[row]; - col_idxs_l[l_diag_idx] = row; - col_idxs_u[u_diag_idx] = row; - vals_l[l_diag_idx] = one(); - vals_u[u_diag_idx] = diag_val; - } + helpers::initialize_l_u( + system_matrix, csr_l, csr_u, + helpers::triangular_mtx_closure([](auto) { return one(); }, + helpers::identity{}), + helpers::triangular_mtx_closure(helpers::identity{}, + helpers::identity{})); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( @@ -248,42 +213,18 @@ void initialize_l(std::shared_ptr exec, const matrix::Csr* system_matrix, matrix::Csr* csr_l, bool diag_sqrt) { - const auto row_ptrs = system_matrix->get_const_row_ptrs(); - const auto col_idxs = system_matrix->get_const_col_idxs(); - const auto vals = system_matrix->get_const_values(); - - const auto row_ptrs_l = csr_l->get_const_row_ptrs(); - auto col_idxs_l = csr_l->get_col_idxs(); - auto vals_l = csr_l->get_values(); - - for (size_type row = 0; row < system_matrix->get_size()[0]; ++row) { - size_type current_index_l = row_ptrs_l[row]; - // if there is no diagonal value, set it to 1 by default - auto diag_val = one(); - for (size_type el = row_ptrs[row]; el < row_ptrs[row + 1]; ++el) { - const auto col = col_idxs[el]; - const auto val = vals[el]; - if (col < row) { - col_idxs_l[current_index_l] = col; - vals_l[current_index_l] = val; - ++current_index_l; - } else if (col == row) { - // save diagonal value - diag_val = val; - } - } - // store diagonal values separately - auto l_diag_idx = row_ptrs_l[row + 1] - 1; - col_idxs_l[l_diag_idx] = row; - // compute square root with sentinel - if (diag_sqrt) { - diag_val = sqrt(diag_val); - if (!is_finite(diag_val)) { - diag_val = one(); - } - } - vals_l[l_diag_idx] = diag_val; - } + helpers::initialize_l(system_matrix, csr_l, + helpers::triangular_mtx_closure( + [diag_sqrt](auto val) { + if (diag_sqrt) { + val = sqrt(val); + if (!is_finite(val)) { + val = one(); + } + } + return val; + }, + helpers::identity{})); } GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( From b59b1a407103a61f01c844201c5acb0761c22b66 Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Tue, 25 Jun 2024 16:07:46 +0200 Subject: [PATCH 2/7] [prec] implement SOR preconditioner reference kernels only reference kernels are available --- common/unified/CMakeLists.txt | 1 + common/unified/preconditioner/sor_kernels.cpp | 44 ++++++ core/CMakeLists.txt | 1 + core/device_hooks/common_kernels.inc.cpp | 11 ++ core/preconditioner/sor.cpp | 133 ++++++++++++++++++ core/preconditioner/sor_kernels.hpp | 50 +++++++ include/ginkgo/core/preconditioner/sor.hpp | 125 ++++++++++++++++ include/ginkgo/ginkgo.hpp | 1 + reference/CMakeLists.txt | 1 + reference/preconditioner/sor_kernels.cpp | 70 +++++++++ 10 files changed, 437 insertions(+) create mode 100644 common/unified/preconditioner/sor_kernels.cpp create mode 100644 core/preconditioner/sor.cpp create mode 100644 core/preconditioner/sor_kernels.hpp create mode 100644 include/ginkgo/core/preconditioner/sor.hpp create mode 100644 reference/preconditioner/sor_kernels.cpp diff --git a/common/unified/CMakeLists.txt b/common/unified/CMakeLists.txt index 00bc21df0c6..132e04c5d9a 100644 --- a/common/unified/CMakeLists.txt +++ b/common/unified/CMakeLists.txt @@ -19,6 +19,7 @@ set(UNIFIED_SOURCES matrix/diagonal_kernels.cpp multigrid/pgm_kernels.cpp preconditioner/jacobi_kernels.cpp + preconditioner/sor_kernels.cpp solver/bicg_kernels.cpp solver/bicgstab_kernels.cpp solver/cg_kernels.cpp diff --git a/common/unified/preconditioner/sor_kernels.cpp b/common/unified/preconditioner/sor_kernels.cpp new file mode 100644 index 00000000000..8932c1df562 --- /dev/null +++ b/common/unified/preconditioner/sor_kernels.cpp @@ -0,0 +1,44 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/preconditioner/sor_kernels.hpp" + +#include +#include + +#include "common/unified/base/kernel_launch.hpp" + + +namespace gko { +namespace kernels { +namespace GKO_DEVICE_NAMESPACE { +namespace sor { + + +template +void initialize_weighted_l( + std::shared_ptr exec, + const matrix::Csr* system_matrix, + remove_complex weight, + matrix::Csr* l_mtx) GKO_NOT_IMPLEMENTED; + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L); + + +template +void initialize_weighted_l_u( + std::shared_ptr exec, + const matrix::Csr* system_matrix, + remove_complex weight, matrix::Csr* l_mtx, + matrix::Csr* u_mtx) GKO_NOT_IMPLEMENTED; + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L_U); + + +} // namespace sor +} // namespace GKO_DEVICE_NAMESPACE +} // namespace kernels +} // namespace gko diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 8c802b2eca5..ef07359e8b4 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -88,6 +88,7 @@ target_sources(${ginkgo_core} multigrid/pgm.cpp multigrid/fixed_coarsening.cpp preconditioner/batch_jacobi.cpp + preconditioner/sor.cpp preconditioner/ic.cpp preconditioner/ilu.cpp preconditioner/isai.cpp diff --git a/core/device_hooks/common_kernels.inc.cpp b/core/device_hooks/common_kernels.inc.cpp index 1ba925e94e3..26de8531741 100644 --- a/core/device_hooks/common_kernels.inc.cpp +++ b/core/device_hooks/common_kernels.inc.cpp @@ -48,6 +48,7 @@ #include "core/preconditioner/batch_jacobi_kernels.hpp" #include "core/preconditioner/isai_kernels.hpp" #include "core/preconditioner/jacobi_kernels.hpp" +#include "core/preconditioner/sor_kernels.hpp" #include "core/reorder/rcm_kernels.hpp" #include "core/solver/batch_bicgstab_kernels.hpp" #include "core/solver/batch_cg_kernels.hpp" @@ -819,6 +820,16 @@ GKO_STUB(GKO_DECLARE_JACOBI_INITIALIZE_PRECISIONS_KERNEL); } // namespace jacobi +namespace sor { + + +GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L); +GKO_STUB_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L_U); + + +} // namespace sor + + namespace isai { diff --git a/core/preconditioner/sor.cpp b/core/preconditioner/sor.cpp new file mode 100644 index 00000000000..30a2539a0cc --- /dev/null +++ b/core/preconditioner/sor.cpp @@ -0,0 +1,133 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include +#include +#include +#include +#include +#include + +#include "core/base/array_access.hpp" +#include "core/base/utils.hpp" +#include "core/factorization/factorization_kernels.hpp" +#include "core/matrix/csr_builder.hpp" +#include "core/preconditioner/sor_kernels.hpp" + +namespace gko { +namespace preconditioner { +namespace { + + +GKO_REGISTER_OPERATION(initialize_row_ptrs_l, + factorization::initialize_row_ptrs_l); +GKO_REGISTER_OPERATION(initialize_row_ptrs_l_u, + factorization::initialize_row_ptrs_l_u); +GKO_REGISTER_OPERATION(initialize_weighted_l, sor::initialize_weighted_l); +GKO_REGISTER_OPERATION(initialize_weighted_l_u, sor::initialize_weighted_l_u); + + +} // namespace + + +template +std::unique_ptr::composition_type> +Sor::generate( + std::shared_ptr system_matrix) const +{ + auto product = + std::unique_ptr(static_cast( + this->LinOpFactory::generate(std::move(system_matrix)).release())); + return product; +} + + +template +std::unique_ptr Sor::generate_impl( + std::shared_ptr system_matrix) const +{ + using Csr = matrix::Csr; + using LTrs = solver::LowerTrs; + using UTrs = solver::UpperTrs; + + GKO_ASSERT_IS_SQUARE_MATRIX(system_matrix); + + auto exec = this->get_executor(); + auto size = system_matrix->get_size(); + + auto csr_matrix = convert_to_with_sorting(exec, system_matrix, + parameters_.skip_sorting); + + 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) { + array l_row_ptrs{exec, size[0] + 1}; + array u_row_ptrs{exec, size[0] + 1}; + exec->run(make_initialize_row_ptrs_l_u( + csr_matrix.get(), l_row_ptrs.get_data(), u_row_ptrs.get_data())); + const auto l_nnz = + static_cast(get_element(l_row_ptrs, size[0])); + const auto u_nnz = + static_cast(get_element(u_row_ptrs, size[0])); + + // create matrices + auto l_mtx = + Csr::create(exec, size, array{exec, l_nnz}, + array{exec, l_nnz}, std::move(l_row_ptrs)); + auto u_mtx = + Csr::create(exec, size, array{exec, u_nnz}, + array{exec, u_nnz}, std::move(u_row_ptrs)); + + // fill l_mtx with 1/w (D + wL) + // fill u_mtx with 1/(1-w) (D + wU) + exec->run(make_initialize_weighted_l_u(csr_matrix.get(), + parameters_.relaxation_factor, + l_mtx.get(), u_mtx.get())); + + // scale u_mtx with D^-1 + auto diag = csr_matrix->extract_diagonal(); + diag->inverse_apply(u_mtx, u_mtx); + + // invert the triangular matrices with triangular solvers + auto l_trs = l_trs_factory->generate(std::move(l_mtx)); + auto u_trs = u_trs_factory->generate(std::move(u_mtx)); + + // return (1/(w * (1 - w)) (D + wL) D^-1 (D + wU))^-1 + // because of the inversion, the factor order is switched + return composition_type::create(std::move(u_trs), std::move(l_trs)); + } else { + array l_row_ptrs{exec, size[0] + 1}; + exec->run(make_initialize_row_ptrs_l(csr_matrix.get(), + l_row_ptrs.get_data())); + const auto l_nnz = + static_cast(get_element(l_row_ptrs, size[0])); + + // create matrices + auto l_mtx = + Csr::create(exec, size, array{exec, l_nnz}, + array{exec, l_nnz}, std::move(l_row_ptrs)); + + // fill l_mtx with 1/w * (D + wL) + exec->run(make_initialize_weighted_l( + csr_matrix.get(), parameters_.relaxation_factor, l_mtx.get())); + + // invert the triangular matrices with triangular solvers + auto l_trs = l_trs_factory->generate(std::move(l_mtx)); + + return composition_type::create(std::move(l_trs)); + } +} + + +#define GKO_DECLARE_SOR(ValueType, IndexType) class Sor + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_SOR); + + +} // namespace preconditioner +} // namespace gko diff --git a/core/preconditioner/sor_kernels.hpp b/core/preconditioner/sor_kernels.hpp new file mode 100644 index 00000000000..fbca3de612c --- /dev/null +++ b/core/preconditioner/sor_kernels.hpp @@ -0,0 +1,50 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_CORE_PRECONDITIONER_SOR_KERNELS_HPP_ +#define GKO_CORE_PRECONDITIONER_SOR_KERNELS_HPP_ + + +#include +#include + +#include "core/base/kernel_declaration.hpp" + + +namespace gko { +namespace kernels { + + +#define GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L(_vtype, _itype) \ + void initialize_weighted_l( \ + std::shared_ptr exec, \ + const matrix::Csr<_vtype, _itype>* system_matrix, \ + remove_complex<_vtype> weight, matrix::Csr<_vtype, _itype>* l_factor) + + +#define GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L_U(_vtype, _itype) \ + void initialize_weighted_l_u( \ + std::shared_ptr exec, \ + const matrix::Csr<_vtype, _itype>* system_matrix, \ + remove_complex<_vtype> weight, matrix::Csr<_vtype, _itype>* l_factor, \ + matrix::Csr<_vtype, _itype>* u_factor) + + +#define GKO_DECLARE_ALL_AS_TEMPLATES \ + template \ + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L(ValueType, IndexType); \ + template \ + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L_U(ValueType, IndexType) + + +GKO_DECLARE_FOR_ALL_EXECUTOR_NAMESPACES(sor, GKO_DECLARE_ALL_AS_TEMPLATES); + + +#undef GKO_DECLARE_ALL_AS_TEMPLATES + + +} // namespace kernels +} // namespace gko + +#endif // GKO_CORE_PRECONDITIONER_SOR_KERNELS_HPP_ diff --git a/include/ginkgo/core/preconditioner/sor.hpp b/include/ginkgo/core/preconditioner/sor.hpp new file mode 100644 index 00000000000..941f012039d --- /dev/null +++ b/include/ginkgo/core/preconditioner/sor.hpp @@ -0,0 +1,125 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_PUBLIC_CORE_PRECONDITIONER_SOR_HPP_ +#define GKO_PUBLIC_CORE_PRECONDITIONER_SOR_HPP_ + + +#include + +#include +#include +#include +#include + + +namespace gko { +namespace preconditioner { + + +/** + * This class generates the (S)SOR preconditioner. + * + * The SOR preconditioner starts from a splitting of the the matrix $A$ into + * $A = D + L + U$, where $L$ contains all entries below the diagonal, and $U$ + * contains all entries above the diagonal. The application of the + * preconditioner is then defined as solving $M x = y$ with + * $$ + * M = \frac{1}{\omega} (D + \omega L), \quad 0 < \omega < 2. + * $$ + * $\omega$ is known as the relaxation factor. + * The preconditioner can be made symmetric, leading to the SSOR preconitioner. + * Here, $M$ is defined as + * $$ + * M = \frac{1}{\omega (2 - \omega)} (D + \omega L) D^{-1} (D + \omega U) , + * \quad 0 < \omega < 2. + * $$ + * + * This class is a factory, which will only generate the preconditioner. The + * resulting LinOp will represent the application of $M^{-1}$. + * + * @tparam ValueType The value type of the internally used CSR matrix + * @tparam IndexType The index type of the internally used CSR matrix + */ +template +class Sor + : public EnablePolymorphicObject, LinOpFactory>, + public EnablePolymorphicAssignment> { + friend class EnablePolymorphicObject; + +public: + struct parameters_type; + friend class enable_parameters_type; + + using value_type = ValueType; + using index_type = IndexType; + using composition_type = Composition; + + struct parameters_type + : public enable_parameters_type { + // skip sorting of input matrix + bool GKO_FACTORY_PARAMETER_SCALAR(skip_sorting, false); + + // determines if SOR or SSOR should be used + bool GKO_FACTORY_PARAMETER_SCALAR(symmetric, false); + + // has to be between 0.0 and 2.0 + remove_complex GKO_FACTORY_PARAMETER_SCALAR( + relaxation_factor, remove_complex(1.2)); + + // factory for the lower triangular factor solver + std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( + l_solver); + + // factory for the upper triangular factor solver, unused if symmetric + // is false + std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( + u_solver); + }; + + /** + * Returns the parameters used to construct the factory. + * + * @return the parameters used to construct the factory. + */ + const parameters_type& get_parameters() { return parameters_; } + + /** + * @copydoc get_parameters + */ + const parameters_type& get_parameters() const { return parameters_; } + + /** + * @copydoc LinOpFactory::generate + * @note This function overrides the default LinOpFactory::generate to + * return a Factorization instead of a generic LinOp, which would need + * to be cast to Factorization again to access its factors. + * It is only necessary because smart pointers aren't covariant. + */ + std::unique_ptr generate( + std::shared_ptr system_matrix) const; + + /** Creates a new parameter_type to set up the factory. */ + static parameters_type build() { return {}; } + +protected: + explicit Sor(std::shared_ptr exec, + const parameters_type& params = {}) + : EnablePolymorphicObject(exec), parameters_(params) + { + GKO_ASSERT(parameters_.relaxation_factor > 0.0 && + parameters_.relaxation_factor < 2.0); + } + + std::unique_ptr generate_impl( + std::shared_ptr system_matrix) const override; + +private: + parameters_type parameters_; +}; +} // namespace preconditioner +} // namespace gko + + +#endif // GKO_PUBLIC_CORE_PRECONDITIONER_SOR_HPP_ diff --git a/include/ginkgo/ginkgo.hpp b/include/ginkgo/ginkgo.hpp index 0fab93dcefe..c44cdee2485 100644 --- a/include/ginkgo/ginkgo.hpp +++ b/include/ginkgo/ginkgo.hpp @@ -118,6 +118,7 @@ #include #include #include +#include #include #include diff --git a/reference/CMakeLists.txt b/reference/CMakeLists.txt index 85b8f33e38b..ab6c210518b 100644 --- a/reference/CMakeLists.txt +++ b/reference/CMakeLists.txt @@ -43,6 +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/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 new file mode 100644 index 00000000000..88ac422dd02 --- /dev/null +++ b/reference/preconditioner/sor_kernels.cpp @@ -0,0 +1,70 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/preconditioner/sor_kernels.hpp" + +#include + +#include +#include + +#include "reference/factorization/factorization_helpers.hpp" + + +namespace gko { +namespace kernels { +namespace reference { +namespace sor { + + +template +void initialize_weighted_l( + std::shared_ptr exec, + const matrix::Csr* system_matrix, + remove_complex weight, matrix::Csr* l_mtx) +{ + auto inv_weight = one(weight) / weight; + factorization::helpers::initialize_l( + system_matrix, l_mtx, + factorization::helpers::triangular_mtx_closure( + [inv_weight](auto val) { return val * inv_weight; }, + [](auto val) { return val; })); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L); + + +template +void initialize_weighted_l_u( + std::shared_ptr exec, + const matrix::Csr* system_matrix, + remove_complex weight, matrix::Csr* l_mtx, + matrix::Csr* u_mtx) +{ + auto inv_weight = one(weight) / weight; + auto inv_two_minus_weight = + one(weight) / (static_cast>(2.0) - weight); + factorization::helpers::initialize_l_u( + system_matrix, l_mtx, u_mtx, + factorization::helpers::triangular_mtx_closure( + [inv_weight](auto val) { return val * inv_weight; }, + factorization::helpers::identity{}), + factorization::helpers::triangular_mtx_closure( + [inv_two_minus_weight](auto val) { + return val * inv_two_minus_weight; + }, + [weight, inv_two_minus_weight](auto val) { + return val * weight * inv_two_minus_weight; + })); +} + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE( + GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L_U); + + +} // namespace sor +} // namespace reference +} // namespace kernels +} // namespace gko From 60be08ab0594541c4bea000aa9f565bd06ec5e41 Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Tue, 25 Jun 2024 17:46:06 +0200 Subject: [PATCH 3/7] [prec] add SOR core + ref tests --- core/test/preconditioner/CMakeLists.txt | 1 + core/test/preconditioner/sor.cpp | 58 +++++ reference/test/preconditioner/CMakeLists.txt | 1 + reference/test/preconditioner/sor_kernels.cpp | 204 ++++++++++++++++++ 4 files changed, 264 insertions(+) create mode 100644 core/test/preconditioner/sor.cpp create mode 100644 reference/test/preconditioner/sor_kernels.cpp diff --git a/core/test/preconditioner/CMakeLists.txt b/core/test/preconditioner/CMakeLists.txt index 41db021e030..87996a79e32 100644 --- a/core/test/preconditioner/CMakeLists.txt +++ b/core/test/preconditioner/CMakeLists.txt @@ -3,3 +3,4 @@ ginkgo_create_test(ic) ginkgo_create_test(ilu) ginkgo_create_test(isai) ginkgo_create_test(jacobi) +ginkgo_create_test(sor) diff --git a/core/test/preconditioner/sor.cpp b/core/test/preconditioner/sor.cpp new file mode 100644 index 00000000000..21c6f1a03e5 --- /dev/null +++ b/core/test/preconditioner/sor.cpp @@ -0,0 +1,58 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include + +#include "core/test/utils.hpp" + + +class SorFactory : public ::testing::Test { +public: + using sor_type = gko::preconditioner::Sor; + using l_isai_type = gko::preconditioner::LowerIsai; + using u_isai_type = gko::preconditioner::UpperIsai; + + std::shared_ptr exec = + gko::ReferenceExecutor::create(); +}; + + +TEST_F(SorFactory, CanDefaultBuild) +{ + auto factory = sor_type::build().on(exec); + + auto params = factory->get_parameters(); + ASSERT_EQ(params.skip_sorting, false); + ASSERT_EQ(params.relaxation_factor, 1.2); + ASSERT_EQ(params.symmetric, false); + ASSERT_EQ(params.l_solver, nullptr); + ASSERT_EQ(params.u_solver, nullptr); +} + + +TEST_F(SorFactory, CanBuildWithParameters) +{ + auto factory = sor_type::build() + .with_skip_sorting(true) + .with_relaxation_factor(0.5) + .with_symmetric(true) + .with_l_solver(l_isai_type::build()) + .with_u_solver(u_isai_type::build()) + .on(exec); + + auto params = factory->get_parameters(); + ASSERT_EQ(params.skip_sorting, true); + ASSERT_EQ(params.relaxation_factor, 0.5); + ASSERT_EQ(params.symmetric, true); + ASSERT_NE(params.l_solver, nullptr); + GKO_ASSERT_DYNAMIC_TYPE(params.l_solver, l_isai_type::Factory); + ASSERT_NE(params.u_solver, nullptr); + GKO_ASSERT_DYNAMIC_TYPE(params.u_solver, u_isai_type::Factory); +} diff --git a/reference/test/preconditioner/CMakeLists.txt b/reference/test/preconditioner/CMakeLists.txt index c8b945e2913..09c88608d65 100644 --- a/reference/test/preconditioner/CMakeLists.txt +++ b/reference/test/preconditioner/CMakeLists.txt @@ -4,3 +4,4 @@ ginkgo_create_test(ic) ginkgo_create_test(isai_kernels) ginkgo_create_test(jacobi) ginkgo_create_test(jacobi_kernels) +ginkgo_create_test(sor_kernels) diff --git a/reference/test/preconditioner/sor_kernels.cpp b/reference/test/preconditioner/sor_kernels.cpp new file mode 100644 index 00000000000..f2bf5f186f9 --- /dev/null +++ b/reference/test/preconditioner/sor_kernels.cpp @@ -0,0 +1,204 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include "core/preconditioner/sor_kernels.hpp" + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" + + +template +class Sor : public ::testing::Test { +public: + using value_type = + typename std::tuple_element<0, decltype(ValueIndexType())>::type; + using index_type = + typename std::tuple_element<1, decltype(ValueIndexType())>::type; + using csr_type = gko::matrix::Csr; + using sor_type = gko::preconditioner::Sor; + + std::shared_ptr exec = + gko::ReferenceExecutor::create(); + gko::remove_complex diag_value = + static_cast>(1.5); + std::shared_ptr mtx = + gko::initialize({{diag_value, 2, 0, 3, 4}, + {-2, diag_value, 5, 0, 0}, + {0, -5, diag_value, 0, 6}, + {-3, 0, 0, diag_value, 7}, + {-4, 0, -6, -7, diag_value}}, + exec); + std::shared_ptr expected_l = + gko::initialize({{diag_value, 0, 0, 0, 0}, + {-2, diag_value, 0, 0, 0}, + {0, -5, diag_value, 0, 0}, + {-3, 0, 0, diag_value, 0}, + {-4, 0, -6, -7, diag_value}}, + exec); + std::shared_ptr expected_u = + gko::initialize({{diag_value, 2, 0, 3, 4}, + {0, diag_value, 5, 0, 0}, + {0, 0, diag_value, 0, 6}, + {0, 0, 0, diag_value, 7}, + {0, 0, 0, 0, diag_value}}, + exec); +}; + +TYPED_TEST_SUITE(Sor, gko::test::ValueIndexTypes, PairTypenameNameGenerator); + + +TYPED_TEST(Sor, CanInitializeLFactor) +{ + using value_type = typename TestFixture::value_type; + auto result = gko::clone(this->expected_l); + result->scale( + gko::initialize>({0.0}, this->exec)); + + gko::kernels::reference::sor::initialize_weighted_l( + this->exec, this->mtx.get(), 1.0, result.get()); + + GKO_ASSERT_MTX_NEAR(result, this->expected_l, 0.0); +} + + +TYPED_TEST(Sor, CanInitializeLFactorWithWeight) +{ + using value_type = typename TestFixture::value_type; + using csr_type = typename TestFixture::csr_type; + auto result = gko::clone(this->expected_l); + 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}}, + this->exec); + + gko::kernels::reference::sor::initialize_weighted_l( + this->exec, this->mtx.get(), this->diag_value, result.get()); + + GKO_ASSERT_MTX_NEAR(result, expected_l, r::value); +} + + +TYPED_TEST(Sor, CanInitializeLAndUFactor) +{ + using value_type = typename TestFixture::value_type; + auto result_l = gko::clone(this->expected_l); + auto result_u = gko::clone(this->expected_u); + result_l->scale( + gko::initialize>({0.0}, this->exec)); + result_u->scale( + gko::initialize>({0.0}, this->exec)); + + gko::kernels::reference::sor::initialize_weighted_l_u( + this->exec, this->mtx.get(), 1.0, result_l.get(), result_u.get()); + + GKO_ASSERT_MTX_NEAR(result_l, this->expected_l, 0.0); + GKO_ASSERT_MTX_NEAR(result_u, this->expected_u, 0.0); +} + + +TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight) +{ + using value_type = typename TestFixture::value_type; + using csr_type = typename TestFixture::csr_type; + auto result_l = gko::clone(this->expected_l); + auto result_u = gko::clone(this->expected_u); + result_l->scale( + 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; + 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}}, + this->exec); + std::shared_ptr expected_u = gko::initialize( + {{this->diag_value * diag_weight, 2 * off_diag_weight, 0, + 3 * off_diag_weight, 4 * off_diag_weight}, + {0, this->diag_value * diag_weight, 5 * off_diag_weight, 0, 0}, + {0, 0, this->diag_value * diag_weight, 0, 6 * off_diag_weight}, + {0, 0, 0, this->diag_value * diag_weight, 7 * off_diag_weight}, + {0, 0, 0, 0, this->diag_value * diag_weight}}, + this->exec); + + gko::kernels::reference::sor::initialize_weighted_l_u( + this->exec, this->mtx.get(), this->diag_value, 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); +} + + +TYPED_TEST(Sor, CanGenerateNonSymmetric) +{ + using value_type = typename TestFixture::value_type; + using index_type = typename TestFixture::index_type; + using sor_type = typename TestFixture::sor_type; + using composition_type = typename sor_type::composition_type; + using trs_type = gko::solver::LowerTrs; + + auto sor_pre = sor_type::build() + .with_relaxation_factor(1.0f) + .on(this->exec) + ->generate(this->mtx); + + testing::StaticAssertTypeEq>(); + const auto& ops = sor_pre->get_operators(); + ASSERT_EQ(ops.size(), 1); + GKO_ASSERT_DYNAMIC_TYPE(ops[0], trs_type); + auto result_l = gko::as(ops[0])->get_system_matrix(); + GKO_ASSERT_MTX_NEAR(result_l, this->expected_l, 0.0); +} + + +TYPED_TEST(Sor, CanGenerateSymmetric) +{ + using value_type = typename TestFixture::value_type; + using index_type = typename TestFixture::index_type; + using sor_type = typename TestFixture::sor_type; + using composition_type = typename sor_type::composition_type; + using l_trs_type = gko::solver::LowerTrs; + using u_trs_type = gko::solver::UpperTrs; + + auto sor_pre = sor_type::build() + .with_symmetric(true) + .with_relaxation_factor(1.0f) + .on(this->exec) + ->generate(this->mtx); + + testing::StaticAssertTypeEq>(); + const auto& ops = sor_pre->get_operators(); + ASSERT_EQ(ops.size(), 2); + GKO_ASSERT_DYNAMIC_TYPE(ops[0], u_trs_type); + GKO_ASSERT_DYNAMIC_TYPE(ops[1], l_trs_type); + auto result_u = gko::as(ops[0])->get_system_matrix(); + auto result_l = gko::as(ops[1])->get_system_matrix(); + GKO_ASSERT_MTX_NEAR(result_l, this->expected_l, 0.0); + auto expected_u = gko::clone(this->expected_u); + expected_u->inv_scale(gko::initialize>( + {this->diag_value}, this->exec)); + GKO_ASSERT_MTX_NEAR(result_u, expected_u, r::value); +} From 47b7166ba59bd072259aea0a10a7f1b583dd9f9e Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Fri, 28 Jun 2024 16:51:24 +0200 Subject: [PATCH 4/7] [prec] add sor parsing --- core/config/config_helper.hpp | 1 + core/config/preconditioner_config.cpp | 2 + core/config/registry.cpp | 1 + core/preconditioner/sor.cpp | 34 +++++++++++++ core/test/config/preconditioner.cpp | 55 +++++++++++++++++++++- include/ginkgo/core/preconditioner/sor.hpp | 6 +++ 6 files changed, 97 insertions(+), 2 deletions(-) diff --git a/core/config/config_helper.hpp b/core/config/config_helper.hpp index 555bb75c2a8..7ddfe35b99a 100644 --- a/core/config/config_helper.hpp +++ b/core/config/config_helper.hpp @@ -64,6 +64,7 @@ enum class LinOpFactoryType : int { Ilu, Isai, Jacobi, + Sor, Multigrid, Pgm, Schwarz diff --git a/core/config/preconditioner_config.cpp b/core/config/preconditioner_config.cpp index cba54cb3356..baf780360e4 100644 --- a/core/config/preconditioner_config.cpp +++ b/core/config/preconditioner_config.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -294,6 +295,7 @@ deferred_factory_parameter parse( GKO_PARSE_VALUE_AND_INDEX_TYPE(Jacobi, gko::preconditioner::Jacobi); +GKO_PARSE_VALUE_AND_INDEX_TYPE(Sor, gko::preconditioner::Sor); } // namespace config diff --git a/core/config/registry.cpp b/core/config/registry.cpp index 188c34b35dd..afba48297ba 100644 --- a/core/config/registry.cpp +++ b/core/config/registry.cpp @@ -44,6 +44,7 @@ configuration_map generate_config_map() {"preconditioner::Ilu", parse}, {"preconditioner::Isai", parse}, {"preconditioner::Jacobi", parse}, + {"preconditioner::Sor", parse}, {"solver::Multigrid", parse}, {"multigrid::Pgm", parse}, #if GINKGO_BUILD_MPI diff --git a/core/preconditioner/sor.cpp b/core/preconditioner/sor.cpp index 30a2539a0cc..0ff534a268c 100644 --- a/core/preconditioner/sor.cpp +++ b/core/preconditioner/sor.cpp @@ -12,6 +12,7 @@ #include "core/base/array_access.hpp" #include "core/base/utils.hpp" +#include "core/config/config_helper.hpp" #include "core/factorization/factorization_kernels.hpp" #include "core/matrix/csr_builder.hpp" #include "core/preconditioner/sor_kernels.hpp" @@ -32,6 +33,39 @@ GKO_REGISTER_OPERATION(initialize_weighted_l_u, sor::initialize_weighted_l_u); } // namespace +template +typename Sor::parameters_type +Sor::parse(const config::pnode& config, + const config::registry& context, + const config::type_descriptor& td_for_child) +{ + auto params = Sor::build(); + + if (auto& obj = config.get("skip_sorting")) { + params.with_skip_sorting(config::get_value(obj)); + } + if (auto& obj = config.get("symmetric")) { + params.with_symmetric(config::get_value(obj)); + } + if (auto& obj = config.get("relaxation_factor")) { + params.with_relaxation_factor( + config::get_value>(obj)); + } + if (auto& obj = config.get("l_solver")) { + params.with_l_solver( + gko::config::parse_or_get_factory( + obj, context, td_for_child)); + } + if (auto& obj = config.get("u_solver")) { + params.with_u_solver( + gko::config::parse_or_get_factory( + obj, context, td_for_child)); + } + + return params; +} + + template std::unique_ptr::composition_type> Sor::generate( diff --git a/core/test/config/preconditioner.cpp b/core/test/config/preconditioner.cpp index 9e81e690967..410e8d74297 100644 --- a/core/test/config/preconditioner.cpp +++ b/core/test/config/preconditioner.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -300,6 +301,55 @@ struct Jacobi }; +struct Sor + : PreconditionerConfigTest<::gko::preconditioner::Sor, + ::gko::preconditioner::Sor> { + using Ir = gko::solver::Ir; + + static pnode::map_type setup_base() + { + return {{"type", pnode{"preconditioner::Sor"}}}; + } + + static void change_template(pnode::map_type& config_map) + { + config_map["value_type"] = pnode{"float32"}; + } + + template + static void set(pnode::map_type& config_map, ParamType& param, registry reg, + std::shared_ptr exec) + { + config_map["skip_sorting"] = pnode{true}; + param.with_skip_sorting(true); + config_map["symmetric"] = pnode{true}; + param.with_symmetric(true); + 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()); + } + + template + static void validate(gko::LinOpFactory* result, AnswerType* answer) + { + auto res_param = gko::as(result)->get_parameters(); + auto ans_param = answer->get_parameters(); + + 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 GINKGO_BUILD_MPI @@ -395,11 +445,12 @@ class Preconditioner : public ::testing::Test { }; -using PreconditionerTypes = ::testing::Types< +using PreconditionerTypes = + ::testing::Types< #if GINKGO_BUILD_MPI ::Schwarz, #endif // GINKGO_BUILD_MPI - ::Ic, ::Ilu, ::Isai, ::Jacobi>; + ::Ic, ::Ilu, ::Isai, ::Jacobi, ::Sor>; TYPED_TEST_SUITE(Preconditioner, PreconditionerTypes, TypenameNameGenerator); diff --git a/include/ginkgo/core/preconditioner/sor.hpp b/include/ginkgo/core/preconditioner/sor.hpp index 941f012039d..276d718dacb 100644 --- a/include/ginkgo/core/preconditioner/sor.hpp +++ b/include/ginkgo/core/preconditioner/sor.hpp @@ -12,6 +12,7 @@ #include #include #include +#include namespace gko { @@ -103,6 +104,11 @@ class Sor /** Creates a new parameter_type to set up the factory. */ static parameters_type build() { return {}; } + static parameters_type parse( + const config::pnode& config, const config::registry& context, + const config::type_descriptor& td_for_child = + config::make_type_descriptor()); + protected: explicit Sor(std::shared_ptr exec, const parameters_type& params = {}) From 23295d024b7d2b79bb9e9ccee384053adb4120ce Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Fri, 28 Jun 2024 14:14:33 +0200 Subject: [PATCH 5/7] =?UTF-8?q?[prec]=20implement=20Gau=C3=9F-Seidel=20pre?= =?UTF-8?q?conditioner?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core/CMakeLists.txt | 3 +- core/preconditioner/gauss_seidel.cpp | 47 +++++++ core/test/preconditioner/CMakeLists.txt | 1 + core/test/preconditioner/gauss_seidel.cpp | 55 ++++++++ .../core/preconditioner/gauss_seidel.hpp | 103 +++++++++++++++ include/ginkgo/ginkgo.hpp | 1 + reference/test/preconditioner/CMakeLists.txt | 1 + .../test/preconditioner/gauss_seidel.cpp | 120 ++++++++++++++++++ 8 files changed, 330 insertions(+), 1 deletion(-) create mode 100644 core/preconditioner/gauss_seidel.cpp create mode 100644 core/test/preconditioner/gauss_seidel.cpp create mode 100644 include/ginkgo/core/preconditioner/gauss_seidel.hpp create mode 100644 reference/test/preconditioner/gauss_seidel.cpp diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index ef07359e8b4..6e0960459e0 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -88,7 +88,8 @@ target_sources(${ginkgo_core} multigrid/pgm.cpp multigrid/fixed_coarsening.cpp preconditioner/batch_jacobi.cpp - preconditioner/sor.cpp + preconditioner/gauss_seidel.cpp + preconditioner/sor.cpp preconditioner/ic.cpp preconditioner/ilu.cpp preconditioner/isai.cpp diff --git a/core/preconditioner/gauss_seidel.cpp b/core/preconditioner/gauss_seidel.cpp new file mode 100644 index 00000000000..20aef490f05 --- /dev/null +++ b/core/preconditioner/gauss_seidel.cpp @@ -0,0 +1,47 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include +#include + + +namespace gko { +namespace preconditioner { + + +template +std::unique_ptr::composition_type> +GaussSeidel::generate( + std::shared_ptr system_matrix) const +{ + auto product = + std::unique_ptr(static_cast( + this->LinOpFactory::generate(std::move(system_matrix)).release())); + return product; +} + + +template +std::unique_ptr GaussSeidel::generate_impl( + std::shared_ptr system_matrix) const +{ + return Sor::build() + .with_skip_sorting(parameters_.skip_sorting) + .with_symmetric(parameters_.symmetric) + .with_relaxation_factor(static_cast>(1.0)) + .with_l_solver(parameters_.l_solver) + .with_u_solver(parameters_.u_solver) + .on(this->get_executor()) + ->generate(std::move(system_matrix)); +} + + +#define GKO_DECLARE_GAUSS_SEIDEL(ValueType, IndexType) \ + class GaussSeidel + +GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GAUSS_SEIDEL); + + +} // namespace preconditioner +} // namespace gko diff --git a/core/test/preconditioner/CMakeLists.txt b/core/test/preconditioner/CMakeLists.txt index 87996a79e32..bb13a60dbe7 100644 --- a/core/test/preconditioner/CMakeLists.txt +++ b/core/test/preconditioner/CMakeLists.txt @@ -1,4 +1,5 @@ ginkgo_create_test(batch_jacobi) +ginkgo_create_test(gauss_seidel) ginkgo_create_test(ic) ginkgo_create_test(ilu) ginkgo_create_test(isai) diff --git a/core/test/preconditioner/gauss_seidel.cpp b/core/test/preconditioner/gauss_seidel.cpp new file mode 100644 index 00000000000..9a2f965db79 --- /dev/null +++ b/core/test/preconditioner/gauss_seidel.cpp @@ -0,0 +1,55 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include + +#include "core/test/utils.hpp" + + +class GaussSeidelFactory : public ::testing::Test { +public: + using GaussSeidel_type = gko::preconditioner::GaussSeidel; + using l_isai_type = gko::preconditioner::LowerIsai; + using u_isai_type = gko::preconditioner::UpperIsai; + + std::shared_ptr exec = + gko::ReferenceExecutor::create(); +}; + + +TEST_F(GaussSeidelFactory, CanDefaultBuild) +{ + auto factory = GaussSeidel_type::build().on(exec); + + auto params = factory->get_parameters(); + ASSERT_EQ(params.skip_sorting, false); + ASSERT_EQ(params.symmetric, false); + ASSERT_EQ(params.l_solver, nullptr); + ASSERT_EQ(params.u_solver, nullptr); +} + + +TEST_F(GaussSeidelFactory, CanBuildWithParameters) +{ + auto factory = GaussSeidel_type::build() + .with_skip_sorting(true) + .with_symmetric(true) + .with_l_solver(l_isai_type::build()) + .with_u_solver(u_isai_type::build()) + .on(exec); + + auto params = factory->get_parameters(); + ASSERT_EQ(params.skip_sorting, true); + ASSERT_EQ(params.symmetric, true); + ASSERT_NE(params.l_solver, nullptr); + GKO_ASSERT_DYNAMIC_TYPE(params.l_solver, l_isai_type::Factory); + ASSERT_NE(params.u_solver, nullptr); + GKO_ASSERT_DYNAMIC_TYPE(params.u_solver, u_isai_type::Factory); +} diff --git a/include/ginkgo/core/preconditioner/gauss_seidel.hpp b/include/ginkgo/core/preconditioner/gauss_seidel.hpp new file mode 100644 index 00000000000..1e482a93819 --- /dev/null +++ b/include/ginkgo/core/preconditioner/gauss_seidel.hpp @@ -0,0 +1,103 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#ifndef GKO_PUBLIC_CORE_PRECONDITIONER_GAUSS_SEIDEL_HPP_ +#define GKO_PUBLIC_CORE_PRECONDITIONER_GAUSS_SEIDEL_HPP_ + + +#include + +#include +#include +#include +#include + + +namespace gko { +namespace preconditioner { + + +/** + * This class generates the Gauss-Seidel preconditioner. + * + * This is the special case of $\omega = 1$ of the (S)SOR preconditioner. + * + * @see Sor + */ +template +class GaussSeidel + : public EnablePolymorphicObject, + LinOpFactory>, + public EnablePolymorphicAssignment> { + friend class EnablePolymorphicObject; + +public: + struct parameters_type; + friend class enable_parameters_type; + + using value_type = ValueType; + using index_type = IndexType; + using composition_type = Composition; + + struct parameters_type + : public enable_parameters_type { + // skip sorting of input matrix + bool GKO_FACTORY_PARAMETER_SCALAR(skip_sorting, false); + + // 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 + std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( + l_solver); + + // factory for the upper triangular factor solver, unused if symmetric + // is false + std::shared_ptr GKO_DEFERRED_FACTORY_PARAMETER( + u_solver); + }; + + /** + * Returns the parameters used to construct the factory. + * + * @return the parameters used to construct the factory. + */ + const parameters_type& get_parameters() { return parameters_; } + + /** + * @copydoc get_parameters + */ + const parameters_type& get_parameters() const { return parameters_; } + + /** + * @copydoc LinOpFactory::generate + * @note This function overrides the default LinOpFactory::generate to + * return a Factorization instead of a generic LinOp, which would need + * to be cast to Factorization again to access its factors. + * It is only necessary because smart pointers aren't covariant. + */ + std::unique_ptr generate( + std::shared_ptr system_matrix) const; + + /** Creates a new parameter_type to set up the factory. */ + static parameters_type build() { return {}; } + +protected: + explicit GaussSeidel(std::shared_ptr exec, + const parameters_type& params = {}) + : EnablePolymorphicObject(exec), + parameters_(params) + {} + + std::unique_ptr generate_impl( + std::shared_ptr system_matrix) const override; + +private: + parameters_type parameters_; +}; +} // namespace preconditioner +} // namespace gko + + +#endif // GKO_PUBLIC_CORE_PRECONDITIONER_GAUSS_SEIDEL_HPP_ diff --git a/include/ginkgo/ginkgo.hpp b/include/ginkgo/ginkgo.hpp index c44cdee2485..61e5b719508 100644 --- a/include/ginkgo/ginkgo.hpp +++ b/include/ginkgo/ginkgo.hpp @@ -114,6 +114,7 @@ #include #include +#include #include #include #include diff --git a/reference/test/preconditioner/CMakeLists.txt b/reference/test/preconditioner/CMakeLists.txt index 09c88608d65..f558aa87495 100644 --- a/reference/test/preconditioner/CMakeLists.txt +++ b/reference/test/preconditioner/CMakeLists.txt @@ -1,4 +1,5 @@ ginkgo_create_test(batch_jacobi_kernels) +ginkgo_create_test(gauss_seidel) ginkgo_create_test(ilu) ginkgo_create_test(ic) ginkgo_create_test(isai_kernels) diff --git a/reference/test/preconditioner/gauss_seidel.cpp b/reference/test/preconditioner/gauss_seidel.cpp new file mode 100644 index 00000000000..2b67b665d77 --- /dev/null +++ b/reference/test/preconditioner/gauss_seidel.cpp @@ -0,0 +1,120 @@ +// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors +// +// SPDX-License-Identifier: BSD-3-Clause + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include "core/test/utils.hpp" +#include "core/utils/matrix_utils.hpp" + + +template +class GaussSeidel : public ::testing::Test { +public: + using value_type = + typename std::tuple_element<0, decltype(ValueIndexType())>::type; + using index_type = + typename std::tuple_element<1, decltype(ValueIndexType())>::type; + using csr_type = gko::matrix::Csr; + using gs_type = gko::preconditioner::GaussSeidel; + using sor_type = gko::preconditioner::Sor; + using ltrs_type = gko::solver::LowerTrs; + using utrs_type = gko::solver::UpperTrs; + + GaussSeidel() + { + auto data = + gko::test::generate_random_matrix_data( + 10, 10, std::uniform_int_distribution<>(2, 6), + std::uniform_real_distribution<>(1, 2), engine); + gko::utils::make_symmetric(data); + gko::utils::make_unit_diagonal(data); + mtx->read(data); + } + + std::default_random_engine engine; + std::shared_ptr exec = + gko::ReferenceExecutor::create(); + std::shared_ptr mtx = csr_type::create(exec); +}; + +TYPED_TEST_SUITE(GaussSeidel, gko::test::ValueIndexTypes, + PairTypenameNameGenerator); + + +TYPED_TEST(GaussSeidel, GenerateSameAsSor) +{ + using real_type = gko::remove_complex; + using gs_type = typename TestFixture::gs_type; + using sor_type = typename TestFixture::sor_type; + using composition_type = typename sor_type::composition_type; + using csr_type = typename TestFixture::csr_type; + using ltrs_type = typename TestFixture::ltrs_type; + + auto gs = gs_type ::build().on(this->exec)->generate(this->mtx); + auto sor = sor_type::build() + .with_relaxation_factor(real_type{1.0}) + .on(this->exec) + ->generate(this->mtx); + + auto gs_comp = dynamic_cast(gs.get()); + auto sor_comp = dynamic_cast(sor.get()); + ASSERT_TRUE(gs_comp); + ASSERT_TRUE(sor_comp); + ASSERT_EQ(gs_comp->get_operators().size(), + sor_comp->get_operators().size()); + GKO_ASSERT_MTX_NEAR( + dynamic_cast(gs_comp->get_operators()[0].get()) + ->get_system_matrix(), + dynamic_cast(sor_comp->get_operators()[0].get()) + ->get_system_matrix(), + 0.0); +} + +TYPED_TEST(GaussSeidel, GenerateSymmetricSameAsSor) +{ + using real_type = gko::remove_complex; + using gs_type = typename TestFixture::gs_type; + using sor_type = typename TestFixture::sor_type; + using composition_type = typename sor_type::composition_type; + using ltrs_type = typename TestFixture::ltrs_type; + using utrs_type = typename TestFixture::utrs_type; + + auto gs = gs_type ::build() + .with_symmetric(true) + .on(this->exec) + ->generate(this->mtx); + auto sor = sor_type::build() + .with_symmetric(true) + .with_relaxation_factor(real_type{1.0}) + .on(this->exec) + ->generate(this->mtx); + + auto gs_comp = dynamic_cast(gs.get()); + auto sor_comp = dynamic_cast(sor.get()); + ASSERT_TRUE(gs_comp); + ASSERT_TRUE(sor_comp); + ASSERT_EQ(gs_comp->get_operators().size(), + sor_comp->get_operators().size()); + GKO_ASSERT_MTX_NEAR( + dynamic_cast(gs_comp->get_operators()[0].get()) + ->get_system_matrix(), + dynamic_cast(sor_comp->get_operators()[0].get()) + ->get_system_matrix(), + 0.0); + GKO_ASSERT_MTX_NEAR( + dynamic_cast(gs_comp->get_operators()[1].get()) + ->get_system_matrix(), + dynamic_cast(sor_comp->get_operators()[1].get()) + ->get_system_matrix(), + 0.0); +} From 17aaedac82f0192112909645516e297e2145d59f Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Tue, 29 Oct 2024 17:51:34 +0000 Subject: [PATCH 6/7] [prec] add gauss seidel parsing --- core/config/config_helper.hpp | 1 + core/config/preconditioner_config.cpp | 2 + core/config/registry.cpp | 2 + core/preconditioner/gauss_seidel.cpp | 31 +++++++++++ core/test/config/preconditioner.cpp | 52 +++++++++++++++++-- .../core/preconditioner/gauss_seidel.hpp | 8 +++ 6 files changed, 93 insertions(+), 3 deletions(-) diff --git a/core/config/config_helper.hpp b/core/config/config_helper.hpp index 7ddfe35b99a..483366765aa 100644 --- a/core/config/config_helper.hpp +++ b/core/config/config_helper.hpp @@ -60,6 +60,7 @@ enum class LinOpFactoryType : int { ParIct, ParIlu, ParIlut, + GaussSeidel, Ic, Ilu, Isai, diff --git a/core/config/preconditioner_config.cpp b/core/config/preconditioner_config.cpp index baf780360e4..68cbf8595ba 100644 --- a/core/config/preconditioner_config.cpp +++ b/core/config/preconditioner_config.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -294,6 +295,7 @@ deferred_factory_parameter parse( } +GKO_PARSE_VALUE_AND_INDEX_TYPE(GaussSeidel, gko::preconditioner::GaussSeidel); GKO_PARSE_VALUE_AND_INDEX_TYPE(Jacobi, gko::preconditioner::Jacobi); GKO_PARSE_VALUE_AND_INDEX_TYPE(Sor, gko::preconditioner::Sor); diff --git a/core/config/registry.cpp b/core/config/registry.cpp index afba48297ba..19da3ed2559 100644 --- a/core/config/registry.cpp +++ b/core/config/registry.cpp @@ -40,6 +40,8 @@ configuration_map generate_config_map() {"factorization::ParIct", parse}, {"factorization::ParIlu", parse}, {"factorization::ParIlut", parse}, + {"preconditioner::GaussSeidel", + parse}, {"preconditioner::Ic", parse}, {"preconditioner::Ilu", parse}, {"preconditioner::Isai", parse}, diff --git a/core/preconditioner/gauss_seidel.cpp b/core/preconditioner/gauss_seidel.cpp index 20aef490f05..aec7a4ff827 100644 --- a/core/preconditioner/gauss_seidel.cpp +++ b/core/preconditioner/gauss_seidel.cpp @@ -5,11 +5,42 @@ #include #include +#include "core/config/config_helper.hpp" + namespace gko { namespace preconditioner { +template +typename GaussSeidel::parameters_type +GaussSeidel::parse( + const config::pnode& config, const config::registry& context, + const config::type_descriptor& td_for_child) +{ + auto params = GaussSeidel::build(); + + if (auto& obj = config.get("skip_sorting")) { + params.with_skip_sorting(config::get_value(obj)); + } + if (auto& obj = config.get("symmetric")) { + params.with_symmetric(config::get_value(obj)); + } + if (auto& obj = config.get("l_solver")) { + params.with_l_solver( + gko::config::parse_or_get_factory( + obj, context, td_for_child)); + } + if (auto& obj = config.get("u_solver")) { + params.with_u_solver( + gko::config::parse_or_get_factory( + obj, context, td_for_child)); + } + + return params; +} + + template std::unique_ptr::composition_type> GaussSeidel::generate( diff --git a/core/test/config/preconditioner.cpp b/core/test/config/preconditioner.cpp index 410e8d74297..c3941f504e2 100644 --- a/core/test/config/preconditioner.cpp +++ b/core/test/config/preconditioner.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -350,6 +351,52 @@ struct Sor }; +struct GaussSeidel + : PreconditionerConfigTest< + ::gko::preconditioner::GaussSeidel, + ::gko::preconditioner::GaussSeidel> { + using Ir = gko::solver::Ir; + + static pnode::map_type setup_base() + { + return {{"type", pnode{"preconditioner::GaussSeidel"}}}; + } + + static void change_template(pnode::map_type& config_map) + { + config_map["value_type"] = pnode{"float32"}; + } + + template + static void set(pnode::map_type& config_map, ParamType& param, registry reg, + std::shared_ptr exec) + { + config_map["skip_sorting"] = pnode{true}; + 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()); + } + + template + static void validate(gko::LinOpFactory* result, AnswerType* answer) + { + auto res_param = gko::as(result)->get_parameters(); + auto ans_param = answer->get_parameters(); + + 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 GINKGO_BUILD_MPI @@ -445,12 +492,11 @@ class Preconditioner : public ::testing::Test { }; -using PreconditionerTypes = - ::testing::Types< +using PreconditionerTypes = ::testing::Types< #if GINKGO_BUILD_MPI ::Schwarz, #endif // GINKGO_BUILD_MPI - ::Ic, ::Ilu, ::Isai, ::Jacobi, ::Sor>; + ::GaussSeidel, ::Ic, ::Ilu, ::Isai, ::Jacobi, ::Sor>; TYPED_TEST_SUITE(Preconditioner, PreconditionerTypes, TypenameNameGenerator); diff --git a/include/ginkgo/core/preconditioner/gauss_seidel.hpp b/include/ginkgo/core/preconditioner/gauss_seidel.hpp index 1e482a93819..7003dd54740 100644 --- a/include/ginkgo/core/preconditioner/gauss_seidel.hpp +++ b/include/ginkgo/core/preconditioner/gauss_seidel.hpp @@ -12,6 +12,7 @@ #include #include #include +#include namespace gko { @@ -83,6 +84,11 @@ class GaussSeidel /** Creates a new parameter_type to set up the factory. */ static parameters_type build() { return {}; } + static parameters_type parse( + const config::pnode& config, const config::registry& context, + const config::type_descriptor& td_for_child = + config::make_type_descriptor()); + protected: explicit GaussSeidel(std::shared_ptr exec, const parameters_type& params = {}) @@ -96,6 +102,8 @@ class GaussSeidel private: parameters_type parameters_; }; + + } // namespace preconditioner } // namespace gko From 4b8387393ff78b9596540195ce75cb0c6df42ad2 Mon Sep 17 00:00:00 2001 From: Marcel Koch Date: Thu, 11 Jul 2024 07:24:09 +0000 Subject: [PATCH 7/7] [sor] review updates: - documentation - tests - don't build upper solver if not symmetric Co-authored-by: Yu-Hsiang M. Tsai --- core/preconditioner/sor.cpp | 5 +-- core/test/config/preconditioner.cpp | 8 ++--- .../core/preconditioner/gauss_seidel.hpp | 7 ++-- include/ginkgo/core/preconditioner/sor.hpp | 6 ++-- reference/CMakeLists.txt | 2 +- reference/test/preconditioner/sor_kernels.cpp | 32 +++++++++---------- 6 files changed, 32 insertions(+), 28 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..c603aaea750 100644 --- a/core/test/config/preconditioner.cpp +++ b/core/test/config/preconditioner.cpp @@ -330,10 +330,10 @@ struct Sor 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()); + param.with_l_solver(Ir::build()); config_map["u_solver"] = pnode{ {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_u_solver(DummyIr::build()); + param.with_u_solver(Ir::build()); } template @@ -377,10 +377,10 @@ struct GaussSeidel param.with_symmetric(true); config_map["l_solver"] = pnode{ {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_l_solver(DummyIr::build()); + param.with_l_solver(Ir::build()); config_map["u_solver"] = pnode{ {{"type", pnode{"solver::Ir"}}, {"value_type", pnode{"float32"}}}}; - param.with_u_solver(DummyIr::build()); + param.with_u_solver(Ir::build()); } template diff --git a/include/ginkgo/core/preconditioner/gauss_seidel.hpp b/include/ginkgo/core/preconditioner/gauss_seidel.hpp index 7003dd54740..75668e652a7 100644 --- a/include/ginkgo/core/preconditioner/gauss_seidel.hpp +++ b/include/ginkgo/core/preconditioner/gauss_seidel.hpp @@ -22,7 +22,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 +50,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); }; diff --git a/include/ginkgo/core/preconditioner/sor.hpp b/include/ginkgo/core/preconditioner/sor.hpp index 276d718dacb..531dded79f2 100644 --- a/include/ginkgo/core/preconditioner/sor.hpp +++ b/include/ginkgo/core/preconditioner/sor.hpp @@ -36,6 +36,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 +71,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); }; 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/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);