[solver] add a saddle point solver (WIP)

dune/xt/la/container/istl.hh

@@ -101,6 +101,8 @@ public:
   using typename ProvidesDataAccess<Traits>::DataType;
   // needed to fix gcc compilation error due to ambiguous lookup of derived type
   using derived_type = typename Traits::derived_type;
+  // for dune-istl's LinearOperator
+  using field_type = ScalarType;
 
 private:
   using MutexesType = typename Traits::MutexesType;

dune/xt/la/solver/istl/amg.hh

@@ -28,47 +28,14 @@
 #include <dune/xt/common/parallel/helper.hh>
 #include <dune/xt/la/container/istl.hh>
 
+#include "preconditioners.hh"
+
 namespace Dune {
 namespace XT {
 namespace LA {
 
 #if HAVE_DUNE_ISTL
 
-template <class O>
-class IdentityPreconditioner : public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
-{
-public:
-  //! \brief The domain type of the preconditioner.
-  typedef typename O::domain_type domain_type;
-  //! \brief The range type of the preconditioner.
-  typedef typename O::range_type range_type;
-  //! \brief The field type of the preconditioner.
-  typedef typename range_type::field_type field_type;
-  typedef O InverseOperator;
-
-  IdentityPreconditioner(const SolverCategory::Category cat)
-    : category_(cat)
-  {}
-
-  //! Category of the preconditioner (see SolverCategory::Category)
-  virtual SolverCategory::Category category() const override final
-  {
-    return category_;
-  }
-
-  virtual void pre(domain_type&, range_type&) override final {}
-
-  virtual void apply(domain_type& v, const range_type& d) override final
-  {
-    v = d;
-  }
-
-  virtual void post(domain_type&) override final {}
-
-private:
-  SolverCategory::Category category_;
-};
-
 //! the general, parallel case
 template <class S, class CommunicatorType>
 class AmgApplicator

dune/xt/la/solver/istl/preconditioners.hh

+// This file is part of the dune-xt-la project:
+//   https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Barbara Verfürth (2015)
+//   Felix Schindler  (2014 - 2017)
+//   Rene Milk        (2014 - 2016, 2018)
+//   Tobias Leibner   (2014, 2017)
+
+#ifndef DUNE_XT_LA_SOLVER_ISTL_PRECONDITIONERS_HH
+#define DUNE_XT_LA_SOLVER_ISTL_PRECONDITIONERS_HH
+
+#include <type_traits>
+#include <cmath>
+
+#if HAVE_DUNE_ISTL
+#  include <dune/istl/preconditioners.hh>
+#endif // HAVE_DUNE_ISTL
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+#if HAVE_DUNE_ISTL
+
+template <class O>
+class IdentityPreconditioner : public Dune::Preconditioner<typename O::domain_type, typename O::range_type>
+{
+public:
+  //! \brief The domain type of the preconditioner.
+  typedef typename O::domain_type domain_type;
+  //! \brief The range type of the preconditioner.
+  typedef typename O::range_type range_type;
+  //! \brief The field type of the preconditioner.
+  typedef typename range_type::field_type field_type;
+  typedef O InverseOperator;
+
+  IdentityPreconditioner(const SolverCategory::Category cat)
+    : category_(cat)
+  {}
+
+  //! Category of the preconditioner (see SolverCategory::Category)
+  virtual SolverCategory::Category category() const override final
+  {
+    return category_;
+  }
+
+  virtual void pre(domain_type&, range_type&) override final {}
+
+  virtual void apply(domain_type& v, const range_type& d) override final
+  {
+    v = d;
+  }
+
+  virtual void post(domain_type&) override final {}
+
+private:
+  SolverCategory::Category category_;
+};
+
+#endif // HAVE_DUNE_ISTL
+
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_SOLVER_ISTL_PRECONDITIONERS_HH

dune/xt/la/solver/istl/saddlepoint.hh

+// This file is part of the dune-xt-la project:
+//   https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Barbara Verfürth (2015)
+//   Felix Schindler  (2014 - 2017)
+//   Rene Milk        (2014 - 2016, 2018)
+//   Tobias Leibner   (2014, 2017)
+
+#ifndef DUNE_XT_LA_SOLVER_ISTL_SADDLEPOINT_HH
+#define DUNE_XT_LA_SOLVER_ISTL_SADDLEPOINT_HH
+
+#include "config.h"
+
+#include <type_traits>
+#include <cmath>
+
+#if HAVE_DUNE_ISTL
+#  include <dune/istl/operators.hh>
+#  include <dune/istl/solvers.hh>
+#endif // HAVE_DUNE_ISTL
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/configuration.hh>
+#include <dune/xt/la/container/istl.hh>
+#include <dune/xt/la/solver.hh>
+
+#include "preconditioners.hh"
+#include "schurcomplement.hh"
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+
+#if HAVE_DUNE_ISTL
+
+
+// Solver for saddle point system (A B1; B2^T C) (u; p) = (f; g) using the Schur complement, i.e., solve (B2^T A^{-1} B1
+// - C) p = B2^T A^{-1} f - g first and then u = A^{-1} (F - B1 p)
+template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+class SaddlePointSolver
+{
+public:
+  using Vector = IstlDenseVector<FieldType>;
+  using Matrix = IstlRowMajorSparseMatrix<FieldType>;
+  using Solver = Solver<Matrix, CommunicatorType>;
+
+  // Matrix and vector dimensions are
+  // A: m x m, B1, B2: m x n, C: n x n, f: m, g: n
+  SaddlePointSolver(const Matrix& A,
+                    const Matrix& B1,
+                    const Matrix& B2,
+                    const Matrix& C,
+                    const Common::Configuration& solver_opts = SolverOptions<Matrix>::options())
+    : schur_complement_(A, B1, B2, C, solver_opts)
+  {}
+
+  static std::vector<std::string> types()
+  {
+    std::vector<std::string> ret{"direct"};
+    ret.insert(ret.begin(), "cg_direct_schurcomplement");
+    return ret;
+  } // ... types()
+
+  static XT::Common::Configuration options(const std::string type = "")
+  {
+    const std::string tp = !type.empty() ? type : types()[0];
+    internal::SolverUtils::check_given(tp, types());
+    Common::Configuration general_opts({"type", "post_check_solves_system", "verbose"}, {tp.c_str(), "1e-5", "0"});
+    Common::Configuration iterative_options({"max_iter", "precision"}, {"10000", "1e-10"});
+    iterative_options += general_opts;
+    if (tp == "direct")
+      return general_opts;
+    else if (tp == "cg_direct_schurcomplement")
+      return iterative_options;
+    else
+      return general_opts;
+  } // ... options(...)
+
+  void apply(const Vector& f, const Vector& g, Vector& u, Vector& p) const
+  {
+    apply(f, g, u, p, types()[0]);
+  }
+
+  void apply(const Vector& f, const Vector& g, Vector& u, Vector& p, const std::string& type) const
+  {
+    apply(f, g, u, p, options(type));
+  }
+
+  int verbosity(const Common::Configuration& opts, const Common::Configuration& default_opts) const
+  {
+    const auto actual_value = opts.get("verbose", default_opts.get<int>("verbose"));
+    return
+#  if HAVE_MPI
+        (communicator_.access().communicator().rank() == 0) ? actual_value : 0;
+#  else
+        actual_value;
+#  endif
+  }
+
+  void apply(const Vector& f, const Vector& g, Vector& u, Vector& p, const Common::Configuration& opts) const
+  {
+    const auto type = opts.get<std::string>("type");
+    if (type == "direct") {
+      // copy matrices to saddle point system matrix
+      // create pattern first
+      const size_t m = f.size();
+      const size_t n = g.size();
+      XT::LA::SparsityPatternDefault system_matrix_pattern(m + n);
+      const auto& A = schur_complement_.A();
+      const auto& B1 = schur_complement_.B1();
+      const auto& B2 = schur_complement_.B2();
+      const auto& C = schur_complement_.C();
+      const auto pattern_A = A.pattern();
+      const auto pattern_B1 = B1.pattern();
+      const auto pattern_B2 = B2.pattern();
+      const auto pattern_C = C.pattern();
+      for (size_t ii = 0; ii < m; ++ii) {
+        for (const auto& jj : pattern_A.inner(ii))
+          system_matrix_pattern.insert(ii, jj);
+        for (const auto& jj : pattern_B1.inner(ii))
+          system_matrix_pattern.insert(ii, m + jj);
+        for (const auto& jj : pattern_B2.inner(ii))
+          system_matrix_pattern.insert(m + jj, ii);
+      } // ii
+      for (size_t ii = 0; ii < n; ++ii)
+        for (const auto& jj : pattern_C.inner(ii))
+          system_matrix_pattern.insert(m + ii, m + jj);
+      system_matrix_pattern.sort();
+
+      // now copy the matrices
+      Matrix system_matrix(m + n, m + n, system_matrix_pattern);
+      for (size_t ii = 0; ii < m; ++ii) {
+        for (const auto& jj : pattern_A.inner(ii))
+          system_matrix.set_entry(ii, jj, A.get_entry(ii, jj));
+        for (const auto& jj : pattern_B1.inner(ii))
+          system_matrix.set_entry(ii, m + jj, B1.get_entry(ii, jj));
+        for (const auto& jj : pattern_B2.inner(ii))
+          system_matrix.set_entry(m + jj, ii, B2.get_entry(ii, jj));
+      } // ii
+      for (size_t ii = 0; ii < n; ++ii)
+        for (const auto& jj : pattern_C.inner(ii))
+          system_matrix.set_entry(m + ii, m + jj, C.get_entry(ii, jj));
+
+      // also copy the rhs
+      Vector system_vector(m + n, 0.), solution_vector(m + n, 0.);
+      for (size_t ii = 0; ii < m; ++ii)
+        system_vector[ii] = f[ii];
+      for (size_t ii = 0; ii < n; ++ii)
+        system_vector[m + ii] = g[ii];
+
+      // solve the system by a direct solver
+      XT::LA::solve(system_matrix, system_vector, solution_vector);
+
+      // copy to result vectors
+      for (size_t ii = 0; ii < m; ++ii)
+        u[ii] = solution_vector[ii];
+      for (size_t ii = 0; ii < n; ++ii)
+        p[ii] = solution_vector[m + ii];
+    } else if (type == "cg_direct_schurcomplement") {
+      // calculate rhs B2^T A^{-1} f - g
+      auto Ainv_f = f;
+      auto rhs_p = g;
+      schur_complement_.A_inv().apply(f, Ainv_f);
+      schur_complement_.B2().mtv(Ainv_f, rhs_p);
+      rhs_p -= g;
+
+      // Solve S p = rhs
+      IdentityPreconditioner<SchurComplementOperator<FieldType, CommunicatorType>> prec(
+          SolverCategory::Category::sequential);
+      auto schur_complement_copy = schur_complement_;
+      Dune::CGSolver<typename Vector::BackendType> outer_solver(schur_complement_copy, prec, 1e-10, 10000, 0, false);
+      InverseOperatorResult res;
+      outer_solver.apply(p.backend(), rhs_p.backend(), res);
+
+      // Now solve u = A^{-1}(f - B1 p)
+      auto rhs_u = f;
+      rhs_u -= schur_complement_.B1() * p;
+      schur_complement_.A_inv().apply(rhs_u, u);
+    }
+  } // ... apply(...)
+
+private:
+  const SchurComplementOperator<FieldType, CommunicatorType> schur_complement_;
+};
+
+#else // HAVE_DUNE_ISTL
+
+template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+class SaddlePointSolver
+{
+  static_assert(Dune::AlwaysFalse<FieldType>::value, "You are missing dune-istl!");
+};
+
+
+#endif // HAVE_DUNE_ISTL
+
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_SOLVER_ISTL_SADDLEPOINT_HH

dune/xt/la/solver/istl/schurcomplement.hh

+// This file is part of the dune-xt-la project:
+//   https://github.com/dune-community/dune-xt-la
+// Copyright 2009-2018 dune-xt-la developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Tobias Leibner   (2019)
+
+#ifndef DUNE_XT_LA_SOLVER_ISTL_SCHURCOMPLEMENT_HH
+#define DUNE_XT_LA_SOLVER_ISTL_SCHURCOMPLEMENT_HH
+
+#if HAVE_DUNE_ISTL
+#  include <dune/istl/operators.hh>
+#  include <dune/istl/solvers.hh>
+#endif // HAVE_DUNE_ISTL
+
+#include <dune/xt/common/exceptions.hh>
+#include <dune/xt/common/configuration.hh>
+#include <dune/xt/la/container/istl.hh>
+#include <dune/xt/la/solver.hh>
+
+namespace Dune {
+namespace XT {
+namespace LA {
+
+
+#if HAVE_DUNE_ISTL
+
+
+// For a saddle point matrix (A B1; B2^T C) this models the Schur complement (B2^T A^{-1} B1 - C)
+template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+class SchurComplementOperator
+  : public Dune::LinearOperator<typename IstlDenseVector<FieldType>::BackendType,
+                                typename IstlDenseVector<FieldType>::BackendType>
+{
+  using BaseType = Dune::LinearOperator<typename IstlDenseVector<FieldType>::BackendType,
+                                        typename IstlDenseVector<FieldType>::BackendType>;
+
+public:
+  using Vector = IstlDenseVector<FieldType>;
+  using VectorBackend = typename Vector::BackendType;
+  using Matrix = IstlRowMajorSparseMatrix<FieldType>;
+  using Solver = Solver<Matrix, CommunicatorType>;
+
+  // Matrix dimensions are
+  // A: m x m, B1, B2: m x n, C: n x n
+  SchurComplementOperator(const Matrix& _A,
+                          const Matrix& _B1,
+                          const Matrix& _B2,
+                          const Matrix& _C,
+                          const Common::Configuration& solver_opts = SolverOptions<Matrix>::options())
+    : A_(_A)
+    , A_inv_(make_solver(A_))
+    , B1_(_B1)
+    , B2_(_B2)
+    , C_(_C)
+    , solver_opts_(solver_opts)
+    , m_vec_1_(_A.rows())
+    , m_vec_2_(_A.rows())
+    , n_vec_1_(_C.rows())
+    , n_vec_2_(_C.rows())
+  {}
+
+  SchurComplementOperator(const SchurComplementOperator& other)
+    : A_(other.A_)
+    , A_inv_(make_solver(A_))
+    , B1_(other.B1_)
+    , B2_(other.B2_)
+    , C_(other.C_)
+    , solver_opts_(other.solver_opts_)
+    , m_vec_1_(other.m_vec_1_)
+    , m_vec_2_(other.m_vec_2_)
+    , n_vec_1_(other.n_vec_1_)
+    , n_vec_2_(other.n_vec_2_)
+  {}
+
+  /*! \brief apply operator to x:  \f$ y = S(x) \f$
+        The input vector is consistent and the output must also be
+     consistent on the interior+border partition.
+   */
+  virtual void apply(const VectorBackend& x, VectorBackend& y) const override final
+  {
+    Vector x_la_vector(x);
+    Vector y_la_vector(y);
+    apply(x_la_vector, y_la_vector);
+    y = y_la_vector.backend();
+  }
+
+  virtual void apply(const Vector& x, Vector& y) const
+  {
+    // we want to calculate y = (B2^T A^{-1} B1 - C) x
+    // calculate B1 x
+    auto& B1x = m_vec_1_;
+    B1_.mv(x, B1x);
+    // calculate A^{-1} B1 x
+    auto& AinvB1x = m_vec_2_;
+    A_inv_.apply(B1x, AinvB1x);
+    // apply B2^T
+    B2_.mtv(AinvB1x, y);
+    // calculate Cx
+    auto& Cx = n_vec_1_;
+    C_.mv(x, Cx);
+    y -= Cx;
+  }
+
+
+  //! apply operator to x, scale and add:  \f$ y = y + \alpha S(x) \f$
+  virtual void applyscaleadd(FieldType alpha, const VectorBackend& x, VectorBackend& y) const override final
+  {
+    Vector x_la_vector(x);
+    Vector y_la_vector(y);
+    applyscaleadd(alpha, x_la_vector, y_la_vector);
+    y = y_la_vector.backend();
+  }
+
+  virtual void applyscaleadd(FieldType alpha, const Vector& x, Vector& y) const
+  {
+    auto Sx = n_vec_2_;
+    apply(x, Sx);
+    Sx *= alpha;
+    y += Sx;
+  }
+
+  //! Category of the linear operator (see SolverCategory::Category)
+  virtual SolverCategory::Category category() const override final
+  {
+    return SolverCategory::Category::sequential;
+  }
+
+  const Solver& A_inv() const
+  {
+    return A_inv_;
+  }
+
+  const Matrix& A() const
+  {
+    return A_;
+  }
+
+  const Matrix& B1() const
+  {
+    return B1_;
+  }
+
+  const Matrix& B2() const
+  {
+    return B2_;
+  }
+
+  const Matrix& C() const
+  {
+    return C_;
+  }
+
+private:
+  const Matrix& A_;
+  const Solver A_inv_;
+  const Matrix& B1_;
+  const Matrix& B2_;
+  const Matrix& C_;
+  const Common::Configuration solver_opts_;
+  // vectors to store intermediate results
+  mutable Vector m_vec_1_;
+  mutable Vector m_vec_2_;
+  mutable Vector n_vec_1_;
+  mutable Vector n_vec_2_;
+};
+
+#else // HAVE_DUNE_ISTL
+
+template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+class SchurComplementOperator
+{
+  static_assert(Dune::AlwaysFalse<FieldType>::value, "You are missing dune-istl!");
+};
+
+#endif // HAVE_DUNE_ISTL
+
+} // namespace LA
+} // namespace XT
+} // namespace Dune
+
+#endif // DUNE_XT_LA_SOLVER_ISTL_SCHURCOMPLEMENT_HH