From 499f37fbd934b6412cbe9dbc70720e4d95b8029e Mon Sep 17 00:00:00 2001
From: Marcel Koch <marcel.koch@kit.edu>
Date: Fri, 28 Jun 2024 14:14:33 +0200
Subject: [PATCH] =?UTF-8?q?[prec]=20implement=20Gau=C3=9F-Seidel=20precond?=
 =?UTF-8?q?itioner?=
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/preconditioner/gauss_seidel.hpp      | 103 ++++++++++++++++++
 include/ginkgo/ginkgo.hpp                     |   1 +
 4 files changed, 153 insertions(+), 1 deletion(-)
 create mode 100644 core/preconditioner/gauss_seidel.cpp
 create mode 100644 include/ginkgo/core/preconditioner/gauss_seidel.hpp

diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt
index aa764e0355d..a8fecaf1bce 100644
--- a/core/CMakeLists.txt
+++ b/core/CMakeLists.txt
@@ -84,7 +84,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 <ginkgo/core/preconditioner/gauss_seidel.hpp>
+#include <ginkgo/core/preconditioner/sor.hpp>
+
+
+namespace gko {
+namespace preconditioner {
+
+
+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
diff --git a/include/ginkgo/core/preconditioner/gauss_seidel.hpp b/include/ginkgo/core/preconditioner/gauss_seidel.hpp
new file mode 100644
index 00000000000..a09e4d93732
--- /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 <vector>
+
+#include <ginkgo/core/base/abstract_factory.hpp>
+#include <ginkgo/core/base/composition.hpp>
+#include <ginkgo/core/base/lin_op.hpp>
+#include <ginkgo/core/base/polymorphic_object.hpp>
+
+
+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 <typename ValueType = default_precision, typename IndexType = int32>
+class GaussSeidel
+    : public EnablePolymorphicObject<GaussSeidel<ValueType, IndexType>,
+                                     LinOpFactory>,
+      public EnablePolymorphicAssignment<GaussSeidel<ValueType, IndexType>> {
+    friend class EnablePolymorphicObject<GaussSeidel, LinOpFactory>;
+
+public:
+    struct parameters_type;
+    friend class enable_parameters_type<parameters_type, GaussSeidel>;
+
+    using value_type = ValueType;
+    using index_type = IndexType;
+    using composition_type = Composition<ValueType>;
+
+    struct parameters_type
+        : public enable_parameters_type<parameters_type, GaussSeidel> {
+        // skip sorting of input matrix
+        bool GKO_FACTORY_PARAMETER_SCALAR(skip_sorting, true);
+
+        // 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<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
+            l_solver);
+
+        // factory for the upper triangular factor solver, unused if symmetric
+        // is false
+        std::shared_ptr<const LinOpFactory> 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<composition_type> generate(
+        std::shared_ptr<const LinOp> system_matrix) const;
+
+    /** Creates a new parameter_type to set up the factory. */
+    static parameters_type build() { return {}; }
+
+protected:
+    explicit GaussSeidel(std::shared_ptr<const Executor> exec,
+                         const parameters_type& params = {})
+        : EnablePolymorphicObject<GaussSeidel, LinOpFactory>(exec),
+          parameters_(params)
+    {}
+
+    std::unique_ptr<LinOp> generate_impl(
+        std::shared_ptr<const LinOp> 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 258891a3712..6078c3305cb 100644
--- a/include/ginkgo/ginkgo.hpp
+++ b/include/ginkgo/ginkgo.hpp
@@ -113,6 +113,7 @@
 #include <ginkgo/core/multigrid/pgm.hpp>
 
 #include <ginkgo/core/preconditioner/batch_jacobi.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>