diff --git a/dune/xt/la/container/container-interface.hh b/dune/xt/la/container/container-interface.hh
index 9c1bf04b2542a469ba1bf62ff0d732be40420d26..a0efbdd58842597c650b18b217ac06f592dad346 100644
--- a/dune/xt/la/container/container-interface.hh
+++ b/dune/xt/la/container/container-interface.hh
@@ -102,6 +102,7 @@ public:
using typename CRTP::derived_type;
using typename CRTP::Traits;
using ScalarType = ScalarImp;
+ using field_type = ScalarType;
using RealType = typename Traits::RealType;
static_assert(std::is_same<ScalarType, typename Traits::ScalarType>::value, "");
diff --git a/dune/xt/la/container/matrix-view.hh b/dune/xt/la/container/matrix-view.hh
index 1ff2ce507f37c24c12b7e05018d0ee04844cfb5b..c8d754b7cd3e30d3da44bffdae27953e1c307ec5 100644
--- a/dune/xt/la/container/matrix-view.hh
+++ b/dune/xt/la/container/matrix-view.hh
@@ -311,6 +311,16 @@ public:
, matrix_(matrix)
{}
+ ThisType& operator=(const Matrix& other)
+ {
+ const auto& patt = const_matrix_view_.get_pattern();
+ assert(patt == other.pattern());
+ for (size_t ii = 0; ii < rows(); ++ii)
+ for (auto&& jj : patt.inner(ii))
+ set_entry(ii, jj, other.get_entry(ii, jj));
+ return *this;
+ }
+
size_t row_index(const size_t ii) const
{
return const_matrix_view_.row_index(ii);
diff --git a/dune/xt/la/container/vector-interface.hh b/dune/xt/la/container/vector-interface.hh
index 27b1ef348d8e22f80212028628e8bd3a8b2f4cfa..48310168dfba40799f5a8d42ed58bd13b68e0a59 100644
--- a/dune/xt/la/container/vector-interface.hh
+++ b/dune/xt/la/container/vector-interface.hh
@@ -303,6 +303,12 @@ public:
return result;
} // ... l1_norm(...)
+ // for compatibility with dune core modules
+ RealType one_norm() const
+ {
+ return l1_norm();
+ }
+
/**
* \brief The l2-norm of the vector.
* \return The l2-norm of the vector.
@@ -316,6 +322,12 @@ public:
// v.dot(v) should always be a ScalarType with zero imaginary part
}
+ // for compatibility with dune core modules
+ RealType two_norm() const
+ {
+ return l2_norm();
+ }
+
/**
* \brief The l-infintiy-norm of the vector.
* \return The l-infintiy-norm of the vector.
@@ -326,6 +338,12 @@ public:
return amax().second;
}
+ // for compatibility with dune core modules
+ RealType inf_norm() const
+ {
+ return sup_norm();
+ }
+
virtual ScalarType standard_deviation() const
{
using std::pow;
diff --git a/dune/xt/la/solver/eigen.hh b/dune/xt/la/solver/eigen.hh
index db67a013f8573e7833121bc63f6c0ad2730d7afd..4c6492cecd901e57687ae4fa92a161460aa50505 100644
--- a/dune/xt/la/solver/eigen.hh
+++ b/dune/xt/la/solver/eigen.hh
@@ -51,21 +51,11 @@ namespace LA {
#if HAVE_EIGEN
-
template <class S, class CommunicatorType>
-class Solver<EigenDenseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
+class SolverOptions<EigenDenseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
{
public:
- typedef EigenDenseMatrix<S> MatrixType;
- typedef typename MatrixType::RealType R;
-
- Solver(const MatrixType& matrix)
- : matrix_(matrix)
- {}
-
- Solver(const MatrixType& matrix, const CommunicatorType& /*communicator*/)
- : matrix_(matrix)
- {}
+ using MatrixType = EigenDenseMatrix<S>;
static std::vector<std::string> types()
{
@@ -76,7 +66,7 @@ public:
"qr.colpivhouseholder",
"qr.fullpivhouseholder",
"lu.fullpiv"};
- } // ... types()
+ } // ... types(...)
static Common::Configuration options(const std::string type = "")
{
@@ -84,11 +74,37 @@ public:
internal::SolverUtils::check_given(tp, types());
Common::Configuration default_options({"type", "post_check_solves_system", "check_for_inf_nan"},
{tp.c_str(), "1e-5", "1"});
- // * for symmetric matrices
+ // for symmetric matrices
if (tp == "ldlt" || tp == "llt") {
default_options.set("pre_check_symmetry", "1e-8");
}
return default_options;
+ }
+};
+
+template <class S, class CommunicatorType>
+class Solver<EigenDenseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
+{
+public:
+ typedef EigenDenseMatrix<S> MatrixType;
+ typedef typename MatrixType::RealType R;
+
+ Solver(const MatrixType& matrix)
+ : matrix_(matrix)
+ {}
+
+ Solver(const MatrixType& matrix, const CommunicatorType& /*communicator*/)
+ : matrix_(matrix)
+ {}
+
+ static std::vector<std::string> types()
+ {
+ return SolverOptions<MatrixType, CommunicatorType>::types();
+ } // ... types()
+
+ static Common::Configuration options(const std::string type = "")
+ {
+ return SolverOptions<MatrixType, CommunicatorType>::options(type);
} // ... options(...)
template <class T1, class T2>
@@ -233,32 +249,11 @@ private:
}; // class Solver
-/**
- * \note lu.sparse will copy the matrix to column major
- * \note qr.sparse will copy the matrix to column major
- * \note ldlt.simplicial will copy the matrix to column major
- * \note llt.simplicial will copy the matrix to column major
- */
template <class S, class CommunicatorType>
-class Solver<EigenRowMajorSparseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
+class SolverOptions<EigenRowMajorSparseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
{
- typedef ::Eigen::SparseMatrix<S, ::Eigen::ColMajor> ColMajorBackendType;
-
public:
- typedef EigenRowMajorSparseMatrix<S> MatrixType;
- typedef typename MatrixType::RealType R;
-
-private:
- typedef typename MatrixType::BackendType::Index EIGEN_size_t;
-
-public:
- Solver(const MatrixType& matrix)
- : matrix_(matrix)
- {}
-
- Solver(const MatrixType& matrix, const CommunicatorType& /*communicator*/)
- : matrix_(matrix)
- {}
+ using MatrixType = EigenRowMajorSparseMatrix<S>;
static std::vector<std::string> types()
{
@@ -291,7 +286,7 @@ public:
// , "superlu" // <- untested
//#endif
};
- } // ... types()
+ } // ... types(...)
static Common::Configuration options(const std::string type = "")
{
@@ -319,6 +314,45 @@ public:
} else if (tp.substr(0, 3) == "cg.")
iterative_options.set("pre_check_symmetry", "1e-8");
return iterative_options;
+ }
+};
+
+
+/**
+ * \note lu.sparse will copy the matrix to column major
+ * \note qr.sparse will copy the matrix to column major
+ * \note ldlt.simplicial will copy the matrix to column major
+ * \note llt.simplicial will copy the matrix to column major
+ */
+template <class S, class CommunicatorType>
+class Solver<EigenRowMajorSparseMatrix<S>, CommunicatorType> : protected internal::SolverUtils
+{
+ typedef ::Eigen::SparseMatrix<S, ::Eigen::ColMajor> ColMajorBackendType;
+
+public:
+ typedef EigenRowMajorSparseMatrix<S> MatrixType;
+ typedef typename MatrixType::RealType R;
+
+private:
+ typedef typename MatrixType::BackendType::Index EIGEN_size_t;
+
+public:
+ Solver(const MatrixType& matrix)
+ : matrix_(matrix)
+ {}
+
+ Solver(const MatrixType& matrix, const CommunicatorType& /*communicator*/)
+ : matrix_(matrix)
+ {}
+
+ static std::vector<std::string> types()
+ {
+ return SolverOptions<MatrixType, CommunicatorType>::types();
+ } // ... types()
+
+ static Common::Configuration options(const std::string type = "")
+ {
+ return SolverOptions<MatrixType, CommunicatorType>::options(type);
} // ... options(...)
template <class T1, class T2>
diff --git a/dune/xt/la/solver/istl/saddlepoint.hh b/dune/xt/la/solver/istl/saddlepoint.hh
index 31a741a9629ef922614959f997f48a08ac1b25e3..348242600941cfb931b1d7ffbab7c1968e868095 100644
--- a/dune/xt/la/solver/istl/saddlepoint.hh
+++ b/dune/xt/la/solver/istl/saddlepoint.hh
@@ -41,12 +41,15 @@ namespace LA {
// 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>
+template <class VectorType = IstlDenseVector<double>,
+ class MatrixType = IstlRowMajorSparseMatrix<double>,
+ class CommunicatorType = SequentialCommunication>
class SaddlePointSolver
{
public:
- using Vector = IstlDenseVector<FieldType>;
- using Matrix = IstlRowMajorSparseMatrix<FieldType>;
+ using Vector = VectorType;
+ using Matrix = MatrixType;
+ using Field = typename VectorType::ScalarType;
// Matrix and vector dimensions are
// A: m x m, B1, B2: m x n, C: n x n, f: m, g: n
@@ -156,7 +159,7 @@ public:
// calculate rhs B2^T A^{-1} f - g
auto Ainv_f = f;
auto rhs_p = g;
- SchurComplementOperator<FieldType, CommunicatorType> schur_complement_op(
+ SchurComplementOperator<Vector, Matrix, CommunicatorType> schur_complement_op(
A_,
B1_,
B2_,
@@ -169,11 +172,11 @@ public:
rhs_p -= g;
// Solve S p = rhs
- IdentityPreconditioner<SchurComplementOperator<FieldType, CommunicatorType>> prec(
+ IdentityPreconditioner<SchurComplementOperator<Vector, Matrix, CommunicatorType>> prec(
SolverCategory::Category::sequential);
- Dune::CGSolver<typename Vector::BackendType> outer_solver(schur_complement_op, prec, 1e-10, 10000, 0, false);
+ Dune::CGSolver<Vector> outer_solver(schur_complement_op, prec, 1e-10, 10000, 0, false);
InverseOperatorResult res;
- outer_solver.apply(p.backend(), rhs_p.backend(), res);
+ outer_solver.apply(p, rhs_p, res);
// Now solve u = A^{-1}(f - B1 p)
auto rhs_u = f;
@@ -191,13 +194,14 @@ private:
#else // HAVE_DUNE_ISTL
-template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+template <class VectorType = IstlDenseVector<double>,
+ class MatrixType = IstlRowMajorSparseMatrix<double>,
+ class CommunicatorType = SequentialCommunication>
class SaddlePointSolver
{
- static_assert(Dune::AlwaysFalse<FieldType>::value, "You are missing dune-istl!");
+ static_assert(Dune::AlwaysFalse<VectorType>::value, "You are missing dune-istl!");
};
-
#endif // HAVE_DUNE_ISTL
} // namespace LA
diff --git a/dune/xt/la/solver/istl/schurcomplement.hh b/dune/xt/la/solver/istl/schurcomplement.hh
index 7a290db42b8c775fe09312aa5df410b5c6202ff7..0f376217e58c22a5172fec79d43c32a07f0a4fb7 100644
--- a/dune/xt/la/solver/istl/schurcomplement.hh
+++ b/dune/xt/la/solver/istl/schurcomplement.hh
@@ -29,18 +29,17 @@ namespace LA {
// 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>
+template <class VectorType = IstlDenseVector<double>,
+ class MatrixType = IstlRowMajorSparseMatrix<double>,
+ class CommunicatorType = SequentialCommunication>
+class SchurComplementOperator : public Dune::LinearOperator<VectorType, VectorType>
{
- using BaseType = Dune::LinearOperator<typename IstlDenseVector<FieldType>::BackendType,
- typename IstlDenseVector<FieldType>::BackendType>;
+ using BaseType = Dune::LinearOperator<VectorType, VectorType>;
public:
- using Vector = IstlDenseVector<FieldType>;
- using VectorBackend = typename Vector::BackendType;
- using Matrix = IstlRowMajorSparseMatrix<FieldType>;
+ using Vector = VectorType;
+ using Matrix = MatrixType;
+ using Field = typename VectorType::ScalarType;
using SolverType = Solver<Matrix, CommunicatorType>;
// Matrix dimensions are
@@ -79,15 +78,7 @@ public:
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
+ virtual void apply(const Vector& x, Vector& y) const override final
{
// we want to calculate y = (B2^T A^{-1} B1 - C) x
// calculate B1 x
@@ -104,17 +95,7 @@ public:
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
+ virtual void applyscaleadd(Field alpha, const Vector& x, Vector& y) const override final
{
auto Sx = n_vec_2_;
apply(x, Sx);
@@ -169,10 +150,12 @@ private:
#else // HAVE_DUNE_ISTL
-template <class FieldType = double, class CommunicatorType = SequentialCommunication>
+template <class VectorType = IstlDenseVector<double>,
+ class MatrixType = IstlRowMajorSparseMatrix<double>,
+ class CommunicatorType = SequentialCommunication>
class SchurComplementOperator
{
- static_assert(Dune::AlwaysFalse<FieldType>::value, "You are missing dune-istl!");
+ static_assert(Dune::AlwaysFalse<VectorType>::value, "You are missing dune-istl!");
};
#endif // HAVE_DUNE_ISTL