Add BlockGS local solver GS and test vs TruncatedBlockGS.

dune/solvers/iterationsteps/blockgssteps.hh

@@ -121,6 +121,29 @@ VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter,
   return x;
 }
 
+static constexpr double defaultGsTol = 9e-11;
+
+//! \brief Solves a block problem using the conjugate gradient method.
+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;
+    }
+  }
+  return x;
+}
+
 } // end namespace LinearSolvers
 
 namespace LocalSolverFromLinearSolver {
@@ -244,6 +267,24 @@ 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.
+template <typename dummy = void>
+auto gs(double tol = LinearSolvers::defaultGsTol,
+        double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
+  return [tol, r](const auto& m, const auto& b, const auto& ignore) {
+    using namespace std::placeholders;
+    auto gsSolver = std::bind(
+        LinearSolvers::gs<std::decay_t<decltype(m)>, std::decay_t<decltype(b)>>,
+        _1, _2, tol);
+    auto gsTruncatedSolver =
+        LocalSolverFromLinearSolver::truncateSymmetrically(gsSolver);
+    auto gsTruncatedRegularizedSolver =
+        LocalSolverRegularizer::diagRegularize(r, gsTruncatedSolver);
+    return gsTruncatedRegularizedSolver(m, b, ignore);
+  };
+}
+
 } // end namespace LocalSolvers
 
 } // end namespace BlockGS

dune/solvers/test/gssteptest.cc

@@ -14,6 +14,7 @@
 #include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
 #include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh>
 #include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 
 #include "common.hh"
 
@@ -108,10 +109,18 @@ struct GSTestSuite {
     }
     {
       BlockGSStep<Matrix, Vector> gsStep1;
-      auto gsStep2 =
-          Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
-              Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
-      diffTest(&gsStep1, "BlockGS", &gsStep2, "Dune::Solvers::DirectBlockGS");
+      auto gsStep2 = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
+          Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+      diffTest(&gsStep1, "BlockGS", &gsStep2, "Dune::Solvers::BlockGS(Direct)");
+    }
+    {
+      TruncatedBlockGSStep<Matrix, Vector> gsStep1(1);
+      auto gsStep2 = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(
+          Dune::Solvers::BlockGS::LocalSolvers::gs(0.0, 0.0));
+      test(&gsStep1, "TruncatedBlockGS recursive");
+      test(&gsStep2, "Dune::Solvers::BlockGS(GS)");
+      diffTest(&gsStep1, "TruncatedBlockGS recursive", &gsStep2,
+               "Dune::Solvers::BlockGS(GS)", 1e-9);
     }
 
     // test regularizable block GS
@@ -132,8 +141,8 @@ struct GSTestSuite {
           SemiDefiniteBlockGSStep<Matrix, Vector> gsStep1(config);
           auto gsString1 = Dune::formatString("SemiDefiniteBlockGS %s %s", type,
                                               reg ? "regularized" : "");
-          auto gsStep2Ptr =
-              Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::createPtrFromConfig(config);
+          auto gsStep2Ptr = Dune::Solvers::BlockGSStepFactory<
+              Matrix, Vector>::createPtrFromConfig(config);
           auto gsString2 = Dune::formatString("Dune::Solvers::BlockGS(%s) %s",
                                               type, reg ? "regularized" : "");
           test(&gsStep1, gsString1);