Add doc and factory for new BlockGS local solver GS.

dune/solvers/iterationsteps/blockgssteps.cc

@@ -33,6 +33,8 @@ namespace Solvers {
       t = BlockGSType::LDLt;
     else if (s == "cg")
       t = BlockGSType::CG;
+    else if (s == "gs")
+      t = BlockGSType::GS;
     else
       lhs.setstate(std::ios_base::failbit);
     return lhs;

dune/solvers/iterationsteps/blockgssteps.hh

@@ -121,25 +121,29 @@ VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter,
   return x;
 }
 
-static constexpr double defaultGsTol = 9e-11;
+static constexpr double defaultGsTol = 0.0;
 
-//! \brief Solves a block problem using the conjugate gradient method.
+/**
+ * \brief Solves a block problem using a single Gauss--Seidel iteration.
+ * \param tol Lines with diagonal entries below this value are skipped and
+ * their respective solutions are set to 0.
+ */
 template <class MBlock, class VBlock>
 VBlock gs(const MBlock& m, const VBlock& b, double tol = defaultGsTol) {
   VBlock x(0);
   for (size_t i = 0; i < m.N(); ++i) {
     const auto& mi = m[i];
     const auto& mii = mi[i];
-    if (std::abs(mii) > tol) {
-      x[i] = b[i];
-      const auto& end = mi.end();
-      for (auto it = mi.begin(); it != end; ++it) {
-        auto j = it.index();
-        if (j != i)
-          x[i] -= (*it) * x[j];
-      }
-      x[i] /= mii;
+    if (std::abs(mii) <= tol)
+      continue;
+    x[i] = b[i];
+    const auto& end = mi.end();
+    for (auto it = mi.begin(); it != end; ++it) {
+      auto j = it.index();
+      if (j != i)
+        x[i] -= (*it) * x[j];
     }
+    x[i] /= mii;
   }
   return x;
 }
@@ -267,8 +271,12 @@ auto cg(size_t maxIter = LinearSolvers::defaultCgMaxIter,
   };
 }
 
-//! \brief is @LinearSolvers::direct with ignore nodes.
-//! \param r Regularization parameter. Set to 0.0 to switch off regularization.
+/**
+ * \brief is @LinearSolvers::gs with ignore nodes.
+ * \param tol Tolerance value for skipping potentially zero diagonal entries
+ * after regularization.
+ * \param r Regularization parameter. Set to 0.0 to switch off regularization.
+ */
 template <typename dummy = void>
 auto gs(double tol = LinearSolvers::defaultGsTol,
         double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
@@ -318,7 +326,7 @@ private:
  * @brief The Type enum provides representatives for (new) block Gauss-Seidel
  *        steps that can be constructed via the @BlockGSStepFactory.
  */
-enum class BlockGSType { Direct, LDLt, CG };
+enum class BlockGSType { Direct, LDLt, CG, GS };
 
 /**
  * @brief Provides support to parse the @BlockGSType enum.
@@ -375,6 +383,11 @@ struct BlockGSStepFactory {
           config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultCgTol);
       return createPtr(BlockGS::LocalSolvers::cg(maxIter, tol, reg), dir);
     }
+    case BlockGSType::GS: {
+      auto tol =
+          config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultGsTol);
+      return createPtr(BlockGS::LocalSolvers::gs(tol, reg), dir);
+    }
     default:
       DUNE_THROW(Dune::NotImplemented,
                  "Factory cannot create this BlockGSType.");