[sor] review updates:

- documentation - tests - don't build upper solver if not symmetric Co-authored-by: Yu-Hsiang M. Tsai <[email protected]>
ginkgo-project · Jul 11, 2024 · db0557d · db0557d
1 parent f8cbc05
commit db0557d
Show file tree

Hide file tree

Showing 6 changed files with 32 additions and 28 deletions.
diff --git a/core/preconditioner/sor.cpp b/core/preconditioner/sor.cpp
@@ -96,10 +96,11 @@ std::unique_ptr<LinOp> Sor<ValueType, IndexType>::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<index_type> l_row_ptrs{exec, size[0] + 1};
         array<index_type> 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
@@ -327,10 +327,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 <bool from_reg, typename AnswerType>
@@ -374,10 +374,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 <bool from_reg, typename AnswerType>

diff --git 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<const LinOpFactory> 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<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
             u_solver);
     };

diff --git 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<value_type> GKO_FACTORY_PARAMETER_SCALAR(
             relaxation_factor, remove_complex<value_type>(1.2));
 
-        // factory for the lower triangular factor solver
+        // factory for the lower triangular factor solver, defaults to LowerTrs
         std::shared_ptr<const LinOpFactory> 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<const LinOpFactory> GKO_DEFERRED_FACTORY_PARAMETER(
             u_solver);
     };

diff --git 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
@@ -80,15 +80,15 @@ TYPED_TEST(Sor, CanInitializeLFactorWithWeight)
     result->scale(
         gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
     std::shared_ptr<csr_type> expected_l =
-        gko::initialize<csr_type>({{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<csr_type>({{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_type>::value);
 }
@@ -122,15 +122,16 @@ TYPED_TEST(Sor, CanInitializeLAndUFactorWithWeight)
         gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
     result_u->scale(
         gko::initialize<gko::matrix::Dense<value_type>>({0.0}, this->exec));
-    auto diag_weight = static_cast<gko::remove_complex<value_type>>(
-        1.0 / (2 - this->diag_value));
-    auto off_diag_weight = this->diag_value * diag_weight;
+    auto factor = static_cast<gko::remove_complex<value_type>>(0.5);
+    auto diag_weight =
+        static_cast<gko::remove_complex<value_type>>(1.0 / (2 - factor));
+    auto off_diag_weight = factor * diag_weight;
     std::shared_ptr<csr_type> expected_l =
-        gko::initialize<csr_type>({{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<csr_type>({{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<csr_type> expected_u = gko::initialize<csr_type>(
         {{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_type>::value);
     GKO_ASSERT_MTX_NEAR(result_u, expected_u, r<value_type>::value);