[prec] implement Gauß-Seidel preconditioner

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

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

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_

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>