Add generalized blockgs steps with modular local solvers.

dune/solvers/iterationsteps/CMakeLists.txt

@@ -2,6 +2,7 @@ install(FILES
     amgstep.hh
     blockgsstep.cc
     blockgsstep.hh
+    blockgssteps.hh
     cgstep.cc
     cgstep.hh
     istlseqilu0step.hh

dune/solvers/iterationsteps/blockgssteps.hh

+#ifndef DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
+#define DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH
+
+#include <algorithm>
+#include <functional>
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/common/staticmatrixtools.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+
+namespace Dune {
+namespace Solvers {
+
+namespace BlockGS {
+
+enum class Direction { FORWARD, BACKWARD, SYMMETRIC };
+
+/**
+ * @brief Provides support to parse the @Direction enum.
+ */
+std::istream& operator>>(std::istream& lhs, Direction& t) {
+  std::string s;
+  lhs >> s;
+  std::transform(s.begin(), s.end(), s.begin(), ::tolower);
+
+  if (s == "forward")
+    t = Direction::FORWARD;
+  else if (s == "backward")
+    t = Direction::BACKWARD;
+  else if (s == "symmetric")
+    t = Direction::SYMMETRIC;
+  else
+    lhs.setstate(std::ios_base::failbit);
+  return lhs;
+}
+
+/**
+ * \brief Iterates over all rows i and updates the i-th component by solving
+ * the i-th correction problem by means of localSolver.
+ * \param m Matrix
+ * \param x Solution vector
+ * \param b Inhomogeneity
+ * \param ignore Ignored DOFs
+ * \param localSolver
+ * \param direction Iteration direction (e.g. forward, backward, ...)
+ */
+template <class M, class V, class BitVector, class LocalSolver>
+void linearStep(const M& m, V& x, const V& b, const BitVector& ignore,
+                LocalSolver&& localSolver,
+                Direction direction = Direction::FORWARD) {
+  // Note: move-capture requires C++14
+  auto blockStep = [&, localSolver = std::move(localSolver) ](size_t i) {
+    const auto& row_i = m[i];
+    // Compute residual
+    auto ri = b[i];
+    for (auto cIt = row_i.begin(); cIt != row_i.end(); ++cIt)
+      cIt->mmv(x[cIt.index()], ri);
+    // Update iterate with correction
+    x[i] += localSolver(row_i[i], std::move(ri), ignore[i]);
+  };
+
+  if (direction != Direction::BACKWARD)
+    for (size_t i = 0; i < x.size(); ++i)
+      blockStep(i);
+  if (direction != Direction::FORWARD)
+    for (size_t i = x.size(); i > 0; --i)
+      blockStep(i - 1);
+}
+
+/// \brief Convenience method for LinearIterationSteps.
+template <class M, class V, class BitVector, class LocalSolver>
+void linearStep(LinearIterationStep<M, V, BitVector>& lis,
+                LocalSolver&& localSolver,
+                Direction direction = Direction::FORWARD) {
+  return linearStep(*lis.mat_, *lis.x_, *lis.rhs_, *lis.ignoreNodes_,
+                    std::forward<LocalSolver>(localSolver), direction);
+}
+
+/**
+ * \brief All functions in this namespace are assumed to compute the solution
+ *        to a correction problem. Hence they do not require an initial guess
+ *        (and assume it to be zero) and return the (approximate) solution.
+ */
+namespace LocalSolvers {
+
+namespace LinearSolvers {
+
+/** \brief Solves a block problem using the MBlock-inherent solve method. See
+ * e.g. FieldMatrix::solve.
+ */
+template <class MBlock, class VBlock>
+VBlock direct(const MBlock& m, const VBlock& b) {
+  VBlock x;
+  m.solve(x, b);
+  return x;
+}
+
+//! \brief Solves a block problem using LDL^t decomposition.
+template <class MBlock, class VBlock>
+VBlock ldlt(MBlock m, const VBlock& b) {
+  VBlock x;
+  StaticMatrix::ldlt(m, m, m);
+  StaticMatrix::ldltSolve(m, m, b, x);
+  return x;
+}
+
+static constexpr size_t defaultCgMaxIter = 100;
+static constexpr double defaultCgTol = 1e-15;
+
+//! \brief Solves a block problem using the conjugate gradient method.
+template <class MBlock, class VBlock>
+VBlock cg(const MBlock& m, VBlock r, size_t maxIter = defaultCgMaxIter,
+          double tol = defaultCgTol) {
+  VBlock x = 0;
+  VBlock p = r;
+  auto q = r;
+  auto rOldSquared = r * r;
+  auto tolSquared = tol * tol;
+
+  for (size_t j = 0; j < maxIter; ++j) {
+    m.mv(p, q);
+    if (rOldSquared < tolSquared)
+      break;
+    auto alpha = rOldSquared / (q * p);
+    x.axpy(alpha, p);
+    r.axpy(-alpha, q);
+    auto rSquared = r * r;
+    auto beta = rSquared / rOldSquared;
+    p *= beta;
+    p += r;
+    rOldSquared = rSquared;
+  }
+  return x;
+}
+
+} // end namespace LinearSolvers
+
+namespace LocalSolverFromLinearSolver {
+// Note: move-capture, auto-return and auto-lambda-arguments require C++14
+/**
+ * \brief Applies truncation to the block problem according to the local ignore
+ * nodes. The system is modified in such a way that for nodes that are to be
+ * ignored, the corresponding rhs is set to 0 while the corresponding matrix row
+ * is set to be the appropriate euclidean basis vector, essentially enforcing an
+ * exact solution for this node to be zero. This leads to zero correction
+ * contributions for the ignored DOFs.
+ * \param ignore The global ignore nodes where the local ignore nodes can be
+ * obtained from.
+ * \param localSolver The block solver used to solve the (modified) local
+ * problem.
+ */
+template <class LinearSolver>
+auto truncate(LinearSolver&& linearSolver) {
+  return [&, localSolver = std::move(linearSolver) ](
+      const auto& m, const auto& b, const auto& ignore) {
+    if (ignore.all())
+      return typename std::result_of<
+          LinearSolver(decltype(m), decltype(b))>::type{0};
+
+    if (ignore.none())
+      return linearSolver(m, b);
+
+    auto mTruncated = m;
+    auto bTruncated = b;
+    assert(b.size() == m.N() && m.N() == m.M());
+    size_t blockSize = b.size();
+    for (size_t j = 0; j < blockSize; ++j) {
+      if (not ignore[j])
+        continue;
+      bTruncated[j] = 0;
+      for (size_t k = 0; k < blockSize; ++k)
+        mTruncated[j][k] = (k == j);
+    }
+    return linearSolver(std::move(mTruncated), std::move(bTruncated));
+  };
+}
+
+} // end namespace LocalSolverFromLinearSolver
+
+namespace LocalSolverRegularizer {
+
+static constexpr double defaultDiagRegularizeParameter = 1e-10;
+
+// Note: move-capture, auto-return and auto-lambda-arguments require C++14
+/**
+ * \brief Applies regularization to the diagonal of the block problem, to
+ *        render possibly indefinite problems definite (and hence enable the
+ *        application of more block problem solve techniques).
+ * \param p Diagonal entries that are non-zero are scaled by (1+p) while
+ *        zero diagonal entries are set to be p.
+ */
+template <class LocalSolver>
+auto diagRegularize(double p, LocalSolver&& localSolver) {
+  return [ p, localSolver = std::move(localSolver) ](
+      const auto& m, const auto& b, const auto& ignore) {
+    if (p == 0)
+      return localSolver(m, b, ignore);
+    auto m_regularized = m;
+    for (size_t j = 0; j < m_regularized.N(); ++j)
+      if (!ignore[j]) {
+        if (m_regularized[j][j] == 0)
+          m_regularized[j][j] = p;
+        else
+          m_regularized[j][j] *= 1.0 + p;
+      }
+    return localSolver(std::move(m_regularized), b, ignore);
+  };
+}
+
+} // end namespace LocalSolverRegularizer
+
+//! \brief is @LinearSolvers::direct with ignore nodes.
+//! \param r Regularization parameter. Set to 0.0 to switch off regularization.
+auto direct(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
+  return [r](const auto& m, const auto& b, const auto& ignore) {
+    auto directSolver = LinearSolvers::direct<std::decay_t<decltype(m)>,
+                                              std::decay_t<decltype(b)>>;
+    auto directTruncatedSolver =
+        LocalSolverFromLinearSolver::truncate(directSolver);
+    auto directTruncatedRegularizedSolver =
+        LocalSolverRegularizer::diagRegularize(r, directTruncatedSolver);
+    return directTruncatedRegularizedSolver(m, b, ignore);
+  };
+}
+
+//! \brief is @LinearSolvers::ldlt with ignore nodes.
+auto ldlt(double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
+  return [r](const auto& m, const auto& b, const auto& ignore) {
+    auto ldltSolver = LinearSolvers::ldlt<std::decay_t<decltype(m)>,
+                                          std::decay_t<decltype(b)>>;
+    auto ldltTruncatedSolver =
+        LocalSolverFromLinearSolver::truncate(ldltSolver);
+    auto ldltTruncatedRegularizedSolver =
+        LocalSolverRegularizer::diagRegularize(r, ldltTruncatedSolver);
+    return ldltTruncatedRegularizedSolver(m, b, ignore);
+  };
+}
+
+//! \brief is @LinearSolvers::cg with ignore nodes.
+auto cg(size_t maxIter = LinearSolvers::defaultCgMaxIter,
+        double tol = LinearSolvers::defaultCgTol,
+        double r = LocalSolverRegularizer::defaultDiagRegularizeParameter) {
+  return [maxIter, tol, r](const auto& m, const auto& b, const auto& ignore) {
+    using namespace std::placeholders;
+    auto cgSolver = std::bind(
+        LinearSolvers::cg<std::decay_t<decltype(m)>, std::decay_t<decltype(b)>>,
+        _1, _2, maxIter, tol);
+    auto cgTruncatedSolver = LocalSolverFromLinearSolver::truncate(cgSolver);
+    auto cgTruncatedRegularizedSolver =
+        LocalSolverRegularizer::diagRegularize(r, cgTruncatedSolver);
+    return cgTruncatedRegularizedSolver(m, b, ignore);
+  };
+}
+
+} // end namespace LocalSolvers
+
+} // end namespace BlockGS
+
+/**
+ * \brief A Gauss--Seidel-type linear iteration step.
+ * \param localSolver The solver how to solve the linear block correction
+ * problems.
+ */
+template <class Matrix, class Vector, class BitVector, class LocalSolver>
+struct BlockGSStep : public LinearIterationStep<Matrix, Vector, BitVector> {
+
+  BlockGSStep(LocalSolver&& localSolver,
+              BlockGS::Direction direction = BlockGS::Direction::FORWARD)
+      : localSolver_(localSolver)
+      , direction_(direction) {}
+
+  void iterate() {
+    assert(this->ignoreNodes_ != nullptr);
+    assert(this->mat_ != nullptr);
+    assert(this->x_ != nullptr);
+    assert(this->rhs_ != nullptr);
+    BlockGS::linearStep(*this, localSolver_, direction_);
+  }
+
+private:
+  LocalSolver localSolver_;
+  BlockGS::Direction direction_;
+};
+
+/**
+ * @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 };
+
+/**
+ * @brief Provides support to parse the @BlockGSType enum.
+ */
+std::istream& operator>>(std::istream& lhs, BlockGSType& t) {
+  std::string s;
+  lhs >> s;
+  std::transform(s.begin(), s.end(), s.begin(), ::tolower);
+
+  if (s == "direct")
+    t = BlockGSType::Direct;
+  else if (s == "ldlt")
+    t = BlockGSType::LDLt;
+  else if (s == "cg")
+    t = BlockGSType::CG;
+  else
+    lhs.setstate(std::ios_base::failbit);
+  return lhs;
+}
+
+/**
+ * @brief The BlockGSStepFactory struct allows to conveniently construct various
+ *        (new) block Gauss-Seidel steps through a ParameterTree config.
+ */
+template <class Matrix, class Vector,
+          class BitVector = DefaultBitVector_t<Vector>>
+struct BlockGSStepFactory {
+  using BlockGSStepPtr =
+      std::shared_ptr<LinearIterationStep<Matrix, Vector, BitVector>>;
+
+  // TODO add ProjectedBlockGSStep (not linear!)
+  // TODO add CachedLDLtBlockGSStep
+
+  //! \brief Create a BlockGSStep with custom local solver.
+  template <class LocalSolver>
+  static auto create(
+      LocalSolver&& localSolver,
+      BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
+    return BlockGSStep<Matrix, Vector, BitVector, LocalSolver>(
+        std::forward<LocalSolver>(localSolver), direction);
+  }
+
+  template <class LocalSolver>
+  static auto createPtr(
+      LocalSolver&& localSolver,
+      BlockGS::Direction direction = BlockGS::Direction::FORWARD) {
+    return std::make_shared<
+        BlockGSStep<Matrix, Vector, BitVector, LocalSolver>>(
+        std::forward<LocalSolver>(localSolver), direction);
+  }
+
+  static BlockGSStepPtr createPtrFromConfig(const Dune::ParameterTree& config) {
+    auto type = config.get<BlockGSType>("type");
+    auto dir = config.get<BlockGS::Direction>("direction",
+                                              BlockGS::Direction::FORWARD);
+    auto reg = config.get("regularize_diag",
+                          BlockGS::LocalSolvers::LocalSolverRegularizer::
+                              defaultDiagRegularizeParameter);
+    switch (type) {
+    case BlockGSType::Direct:
+      return createPtr(BlockGS::LocalSolvers::direct(reg), dir);
+    case BlockGSType::LDLt:
+      return createPtr(BlockGS::LocalSolvers::ldlt(reg), dir);
+    case BlockGSType::CG: {
+      auto maxIter = config.get(
+          "maxIter", BlockGS::LocalSolvers::LinearSolvers::defaultCgMaxIter);
+      auto tol =
+          config.get("tol", BlockGS::LocalSolvers::LinearSolvers::defaultCgTol);
+      return createPtr(BlockGS::LocalSolvers::cg(maxIter, tol, reg), dir);
+    }
+    default:
+      DUNE_THROW(Dune::NotImplemented,
+                 "Factory cannot create this BlockGSType.");
+    }
+  }
+};
+
+} // end namespace Solvers
+} // end namespace Dune
+
+#endif // DUNE_SOLVERS_ITERATIONSTEPS_BLOCKGSSTEPS_HH

dune/solvers/test/gssteptest.cc

@@ -13,6 +13,7 @@
 #include <dune/solvers/iterationsteps/blockgsstep.hh>
 #include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
 #include <dune/solvers/iterationsteps/semidefiniteblockgsstep.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
 
 #include "common.hh"
 
@@ -22,11 +23,9 @@
   *  is solved correctly for a random rhs by a LoopSolver employing
   *  a GSStep.
   */
-struct GSTestSuite
-{
+struct GSTestSuite {
   template <class GridType>
-  bool check(const GridType& grid)
-  {
+  bool check(const GridType& grid) {
     bool passed = true;
 
     using Problem = SymmetricSampleProblem<1, typename GridType::LevelGridView>;
@@ -37,6 +36,17 @@ struct GSTestSuite
 
     using Matrix = typename Problem::Matrix;
     using Vector = typename Problem::Vector;
+    enum { BlockSize = Vector::block_type::dimension };
+
+    // create some obstacle
+    size_t size = p.u.size();
+    std::vector<BoxConstraint<double, BlockSize>> obstacles(size);
+    for (size_t i = 0; i < size; ++i) {
+      for (size_t j = 0; j < BlockSize; ++j) {
+        obstacles[i].lower(j) = j;
+        obstacles[i].upper(j) = j + 1;
+      }
+    }
 
     using LoopSolver = ::LoopSolver<Vector>;
     using Step = LinearIterationStep<Matrix, Vector>;
@@ -46,8 +56,8 @@ struct GSTestSuite
       step->setProblem(p.A, u_copy, p.rhs);
       step->ignoreNodes_ = &p.ignore;
 
-      LoopSolver solver(step, maxIterations, stepTol,
-                        &p.energyNorm, verbosity, relativeErrors);
+      LoopSolver solver(step, maxIterations, stepTol, &p.energyNorm, verbosity,
+                        relativeErrors);
       solver.check();
       solver.preprocess();
       solver.solve();
@@ -55,11 +65,12 @@ struct GSTestSuite
       return u_copy;
     };
 
-    auto analyse = [&](const Vector& v, std::string testCase, double tolerance) {
+    auto analyse =
+        [&](const Vector& v, std::string testCase, double tolerance) {
       const auto normDiff = p.energyNorm.diff(p.u_ex, v);
       const std::string s = Dune::formatString("%*s", -60, testCase.c_str());
-      std::cout << "error for solver \"" << s
-                << "\": " << normDiff << std::endl;
+      std::cout << "error for solver \"" << s << "\": " << normDiff
+                << std::endl;
       if (normDiff > tolerance) {
         std::cerr << "Error too large for solver \"" << testCase << "\"."
                   << std::endl;
@@ -73,25 +84,63 @@ struct GSTestSuite
       passed = passed and analyse(result, name, 1e-7);
     };
 
+    auto diffTest = [&](Step* step1, std::string name1, Step* step2,
+                        std::string name2, double maxDiff = 0.0) {
+      auto result1 = solve(step1, 0.0, 5);
+      auto result2 = solve(step2, 0.0, 5);
+
+      result1 -= result2;
+      if (result1.two_norm() > maxDiff) {
+        passed = false;
+        std::cerr << name1 << " and " << name2 << " differ with norm "
+                  << result1.two_norm() << "." << std::endl;
+        return;
+      }
+      std::cout << Dune::formatString("%*s", -60, name1.c_str()) << " and "
+                << Dune::formatString("%*s", -60, name2.c_str()) << " coincide"
+                << ((maxDiff > 0.0)
+                        ? " up to " + Dune::formatString("%e", maxDiff)
+                        : ".") << std::endl;
+      passed = passed and true;
+    };
+
     {
       BlockGSStep<Matrix, Vector> gsStep;
       test(&gsStep, "BlockGS");
     }
+    {
+      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");
+    }
 
     // test regularizable block GS
-    for(auto reg : {false, true}) {
-      for(auto type : {"Direct", "LDLt", "CG"}) {
+    for (auto reg : {false, true}) {
+      for (auto type : {"Direct", "LDLt", "CG"}) {
         Dune::ParameterTree config;
         config["type"] = type;
-        config["regularize_diagonal"] = reg ? "1" : "0"; // for SemiDefiniteBlockGS
-        config["regularize_diag"] = reg ? "1e-10" : "0.0"; // for Dune::Solvers::*BlockGS
-        config["maxIter"] = "100"; // for Dune::Solvers::CGBlockGS
-        config["tol"] = "1e-10"; // for Dune::Solvers::CGBlockGS
+        config["regularize_diagonal"] =
+            reg ? "1" : "0"; // for SemiDefiniteBlockGS
+        config["regularize_diag"] =
+            reg ? "1e-10" : "0.0";    // for Dune::Solvers::*BlockGS
+        config["cg_maxIter"] = "100"; // for SemiDefiniteBlockGS with CG
+        config["cg_tol"] = "1e-10";   // for Dune::Solvers::CGBlockGS with CG
+        config["maxIter"] = "100";    // for Dune::Solvers::CGBlockGS
+        config["tol"] = "1e-10";      // for Dune::Solvers::CGBlockGS
 
         {
-          SemiDefiniteBlockGSStep<Matrix, Vector> gsStep(config);
-          test(&gsStep, Dune::formatString("SemiDefiniteBlockGS %s %s",
-                                           type, reg ? "regularized" : ""));
+          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 gsString2 = Dune::formatString("Dune::Solvers::BlockGS(%s) %s",
+                                              type, reg ? "regularized" : "");
+          test(&gsStep1, gsString1);
+          test(gsStep2Ptr.get(), gsString2);
+          diffTest(&gsStep1, gsString1, gsStep2Ptr.get(), gsString2, 1e-7);
         }
       }
     }
@@ -107,14 +156,6 @@ struct GSTestSuite
         test(&gsStep, "ProjectedBlockGS free");
       }
       hasObstacle.setAll();
-      std::vector<typename GSStep::Obstacle> obstacles(size);
-      for(size_t i=0; i<size; ++i) {
-        hasObstacle[i] = true;
-        for(size_t j=0; j<Vector::block_type::dimension; ++j) {
-          obstacles[i].lower(j) = j;
-          obstacles[i].upper(j) = j+1;
-        }
-      }
       {
         ProjectedBlockGSStep<Matrix, Vector> gsStep;
         gsStep.hasObstacle_ = &hasObstacle;
@@ -127,8 +168,7 @@ struct GSTestSuite
   }
 };
 
-int main(int argc, char** argv)
-{
+int main(int argc, char** argv) {
   Dune::MPIHelper::instance(argc, argv);
   bool passed(true);