Adds (S)SOR Preconditioner

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

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 <ginkgo/core/base/math.hpp>
+#include <ginkgo/core/matrix/csr.hpp>
+
+#include "common/unified/base/kernel_launch.hpp"
+
+
+namespace gko {
+namespace kernels {
+namespace GKO_DEVICE_NAMESPACE {
+namespace sor {
+
+
+template <typename ValueType, typename IndexType>
+void initialize_weighted_l(
+    std::shared_ptr<const DefaultExecutor> exec,
+    const matrix::Csr<ValueType, IndexType>* system_matrix,
+    remove_complex<ValueType> weight,
+    matrix::Csr<ValueType, IndexType>* l_mtx) GKO_NOT_IMPLEMENTED;
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(
+    GKO_DECLARE_SOR_INITIALIZE_WEIGHTED_L);
+
+
+template <typename ValueType, typename IndexType>
+void initialize_weighted_l_u(
+    std::shared_ptr<const DefaultExecutor> exec,
+    const matrix::Csr<ValueType, IndexType>* system_matrix,
+    remove_complex<ValueType> weight, matrix::Csr<ValueType, IndexType>* l_mtx,
+    matrix::Csr<ValueType, IndexType>* 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

core/CMakeLists.txt

@@ -88,6 +88,8 @@ target_sources(${ginkgo_core}
     multigrid/pgm.cpp
     multigrid/fixed_coarsening.cpp
     preconditioner/batch_jacobi.cpp
+    preconditioner/gauss_seidel.cpp
+    preconditioner/sor.cpp
     preconditioner/ic.cpp
     preconditioner/ilu.cpp
     preconditioner/isai.cpp

core/config/config_helper.hpp

@@ -60,10 +60,12 @@ enum class LinOpFactoryType : int {
     ParIct,
     ParIlu,
     ParIlut,
+    GaussSeidel,
     Ic,
     Ilu,
     Isai,
     Jacobi,
+    Sor,
     Multigrid,
     Pgm,
     Schwarz

core/config/preconditioner_config.cpp

@@ -5,10 +5,12 @@
 #include <ginkgo/core/base/exception_helpers.hpp>
 #include <ginkgo/core/config/config.hpp>
 #include <ginkgo/core/config/registry.hpp>
+#include <ginkgo/core/preconditioner/gauss_seidel.hpp>
 #include <ginkgo/core/preconditioner/ic.hpp>
 #include <ginkgo/core/preconditioner/ilu.hpp>
 #include <ginkgo/core/preconditioner/isai.hpp>
 #include <ginkgo/core/preconditioner/jacobi.hpp>
+#include <ginkgo/core/preconditioner/sor.hpp>
 #include <ginkgo/core/solver/gmres.hpp>
 #include <ginkgo/core/solver/ir.hpp>
 #include <ginkgo/core/solver/triangular.hpp>
@@ -293,7 +295,9 @@ deferred_factory_parameter<gko::LinOpFactory> parse<LinOpFactoryType::Isai>(
 }
 
 
+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);
 
 
 }  // namespace config

core/config/registry.cpp

@@ -40,10 +40,13 @@ configuration_map generate_config_map()
             {"factorization::ParIct", parse<LinOpFactoryType::ParIct>},
             {"factorization::ParIlu", parse<LinOpFactoryType::ParIlu>},
             {"factorization::ParIlut", parse<LinOpFactoryType::ParIlut>},
+            {"preconditioner::GaussSeidel",
+             parse<LinOpFactoryType::GaussSeidel>},
             {"preconditioner::Ic", parse<LinOpFactoryType::Ic>},
             {"preconditioner::Ilu", parse<LinOpFactoryType::Ilu>},
             {"preconditioner::Isai", parse<LinOpFactoryType::Isai>},
             {"preconditioner::Jacobi", parse<LinOpFactoryType::Jacobi>},
+            {"preconditioner::Sor", parse<LinOpFactoryType::Sor>},
             {"solver::Multigrid", parse<LinOpFactoryType::Multigrid>},
             {"multigrid::Pgm", parse<LinOpFactoryType::Pgm>},
 #if GINKGO_BUILD_MPI

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 @@
 }  // 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 {
 
 

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 <utility>
+
+
+namespace gko {
+namespace factorization {
+
+
+struct identity {
+    template <typename T>
+    constexpr T operator()(T value)
+    {
+        return value;
+    }
+};
+
+
+template <typename DiagClosure, typename OffDiagClosure>
+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 <typename T>
+    constexpr T map_diag(T value)
+    {
+        return diag_closure_(value);
+    }
+
+    template <typename T>
+    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

core/preconditioner/gauss_seidel.cpp

@@ -0,0 +1,78 @@
+// SPDX-FileCopyrightText: 2017 - 2024 The Ginkgo authors
+//
+// SPDX-License-Identifier: BSD-3-Clause
+
+#include <ginkgo/core/preconditioner/gauss_seidel.hpp>
+#include <ginkgo/core/preconditioner/sor.hpp>
+
+#include "core/config/config_helper.hpp"
+
+
+namespace gko {
+namespace preconditioner {
+
+
+template <typename ValueType, typename IndexType>
+typename GaussSeidel<ValueType, IndexType>::parameters_type
+GaussSeidel<ValueType, IndexType>::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<bool>(obj));
+    }
+    if (auto& obj = config.get("symmetric")) {
+        params.with_symmetric(config::get_value<bool>(obj));
+    }
+    if (auto& obj = config.get("l_solver")) {
+        params.with_l_solver(
+            gko::config::parse_or_get_factory<const LinOpFactory>(
+                obj, context, td_for_child));
+    }
+    if (auto& obj = config.get("u_solver")) {
+        params.with_u_solver(
+            gko::config::parse_or_get_factory<const LinOpFactory>(
+                obj, context, td_for_child));
+    }
+
+    return params;
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<typename GaussSeidel<ValueType, IndexType>::composition_type>
+GaussSeidel<ValueType, IndexType>::generate(
+    std::shared_ptr<const LinOp> system_matrix) const
+{
+    auto product =
+        std::unique_ptr<composition_type>(static_cast<composition_type*>(
+            this->LinOpFactory::generate(std::move(system_matrix)).release()));
+    return product;
+}
+
+
+template <typename ValueType, typename IndexType>
+std::unique_ptr<LinOp> GaussSeidel<ValueType, IndexType>::generate_impl(
+    std::shared_ptr<const LinOp> system_matrix) const
+{
+    return Sor<ValueType, IndexType>::build()
+        .with_skip_sorting(parameters_.skip_sorting)
+        .with_symmetric(parameters_.symmetric)
+        .with_relaxation_factor(static_cast<remove_complex<ValueType>>(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<ValueType, IndexType>
+
+GKO_INSTANTIATE_FOR_EACH_VALUE_AND_INDEX_TYPE(GKO_DECLARE_GAUSS_SEIDEL);
+
+
+}  // namespace preconditioner
+}  // namespace gko