From 3fda7d9a81a44708fa7ba928c3933bb8ae26c69c Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Fri, 1 Feb 2019 15:04:50 +0100
Subject: [PATCH] [solver] add a saddle point solver (WIP)
---
dune/xt/la/container/istl.hh | 2 +
dune/xt/la/solver/istl/amg.hh | 37 +---
dune/xt/la/solver/istl/preconditioners.hh | 70 ++++++++
dune/xt/la/solver/istl/saddlepoint.hh | 206 ++++++++++++++++++++++
dune/xt/la/solver/istl/schurcomplement.hh | 184 +++++++++++++++++++
5 files changed, 464 insertions(+), 35 deletions(-)
create mode 100644 dune/xt/la/solver/istl/preconditioners.hh
create mode 100644 dune/xt/la/solver/istl/saddlepoint.hh
create mode 100644 dune/xt/la/solver/istl/schurcomplement.hh
diff --git a/dune/xt/la/container/istl.hh b/dune/xt/la/container/istl.hh
index ed7f4f2a0..884671707 100644
--- a/dune/xt/la/container/istl.hh
+++ b/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;
diff --git a/dune/xt/la/solver/istl/amg.hh b/dune/xt/la/solver/istl/amg.hh
index fe81068cd..ceded1782 100644
--- a/dune/xt/la/solver/istl/amg.hh
+++ b/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
diff --git a/dune/xt/la/solver/istl/preconditioners.hh b/dune/xt/la/solver/istl/preconditioners.hh
new file mode 100644
index 000000000..f35112d4d
--- /dev/null
+++ b/dune/xt/la/solver/istl/preconditioners.hh
@@ -0,0 +1,70 @@
+// 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
diff --git a/dune/xt/la/solver/istl/saddlepoint.hh b/dune/xt/la/solver/istl/saddlepoint.hh
new file mode 100644
index 000000000..8f76bf3b8
--- /dev/null
+++ b/dune/xt/la/solver/istl/saddlepoint.hh
@@ -0,0 +1,206 @@
+// 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
diff --git a/dune/xt/la/solver/istl/schurcomplement.hh b/dune/xt/la/solver/istl/schurcomplement.hh
new file mode 100644
index 000000000..4bdd234c4
--- /dev/null
+++ b/dune/xt/la/solver/istl/schurcomplement.hh
@@ -0,0 +1,184 @@
+// 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
--
GitLab