diff --git a/dune/xt/la/container/common.hh b/dune/xt/la/container/common.hh
index 8ae76d14a15e9e31f9b4a51d97dde0ac9067d9c8..01dbb6ede6c2d0912dc2655a33a03c671d7988f6 100644
--- a/dune/xt/la/container/common.hh
+++ b/dune/xt/la/container/common.hh
@@ -632,6 +632,11 @@ public:
/// \}
+ using MatrixInterfaceType::operator+;
+ using MatrixInterfaceType::operator-;
+ using MatrixInterfaceType::operator+=;
+ using MatrixInterfaceType::operator-=;
+
protected:
/**
* \see ContainerInterface
@@ -944,6 +949,18 @@ public:
/// \}
+ using MatrixInterfaceType::operator+;
+ using MatrixInterfaceType::operator-;
+ using MatrixInterfaceType::operator+=;
+ using MatrixInterfaceType::operator-=;
+
+protected:
+ inline void ensure_uniqueness()
+ {
+ if (!entries_.unique())
+ entries_ = std::make_shared<EntriesVectorType>(*entries_);
+ } // ... ensure_uniqueness(...)
+
private:
size_t get_entry_index(const size_t rr, const size_t cc, const bool throw_if_not_in_pattern = true) const
{
@@ -959,14 +976,6 @@ private:
return size_t(-1);
}
-protected:
- inline void ensure_uniqueness()
- {
- if (!entries_.unique())
- entries_ = std::make_shared<EntriesVectorType>(*entries_);
- } // ... ensure_uniqueness(...)
-
-private:
size_t num_rows_, num_cols_;
std::shared_ptr<EntriesVectorType> entries_;
std::shared_ptr<IndexVectorType> row_pointers_;
diff --git a/dune/xt/la/container/eigen/dense.hh b/dune/xt/la/container/eigen/dense.hh
index 5f261be8f732a43f20450ef49fb4fd93faae0023..cf926cff3a67c07281514b5a18766d3cffe5e462 100644
--- a/dune/xt/la/container/eigen/dense.hh
+++ b/dune/xt/la/container/eigen/dense.hh
@@ -647,6 +647,11 @@ public:
* \}
*/
+ using MatrixInterfaceType::operator+;
+ using MatrixInterfaceType::operator-;
+ using MatrixInterfaceType::operator+=;
+ using MatrixInterfaceType::operator-=;
+
protected:
inline void ensure_uniqueness()
{
diff --git a/dune/xt/la/container/eigen/sparse.hh b/dune/xt/la/container/eigen/sparse.hh
index 281f7fbbb98e8b7a9d7bcbf94fe24f2ab3f13a3e..799b42bb5dced8e0325568385043e057e3afc397 100644
--- a/dune/xt/la/container/eigen/sparse.hh
+++ b/dune/xt/la/container/eigen/sparse.hh
@@ -434,6 +434,11 @@ public:
/// \}
+ using MatrixInterfaceType::operator+;
+ using MatrixInterfaceType::operator-;
+ using MatrixInterfaceType::operator+=;
+ using MatrixInterfaceType::operator-=;
+
private:
bool these_are_valid_indices(const size_t ii, const size_t jj) const
{
diff --git a/dune/xt/la/container/istl.hh b/dune/xt/la/container/istl.hh
index fae0045618124f4f13e88013803aff2667bf358f..fac0a369c96a617c2c3e644a9fa63a38671b74f1 100644
--- a/dune/xt/la/container/istl.hh
+++ b/dune/xt/la/container/istl.hh
@@ -690,6 +690,11 @@ public:
/// \}
+ using MatrixInterfaceType::operator+;
+ using MatrixInterfaceType::operator-;
+ using MatrixInterfaceType::operator+=;
+ using MatrixInterfaceType::operator-=;
+
private:
void build_sparse_matrix(const size_t rr, const size_t cc, const SparsityPatternDefault& patt)
{
diff --git a/dune/xt/la/container/matrix-interface.hh b/dune/xt/la/container/matrix-interface.hh
index a3c4b6bcc1effaca089131706675f85d785cf8fb..62c1037cf88c5b71a735d025b4d0e6c5cd893075 100644
--- a/dune/xt/la/container/matrix-interface.hh
+++ b/dune/xt/la/container/matrix-interface.hh
@@ -114,6 +114,36 @@ public:
CHECK_AND_CALL_CRTP(this->as_imp().unit_col(jj));
}
+ template <class MM>
+ derived_type operator*(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ return multiply(other);
+ }
+
+ template <class MM>
+ derived_type operator+(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ return add(other);
+ }
+
+ template <class MM>
+ derived_type operator-(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ return subtract(other);
+ }
+
+ template <class MM>
+ derived_type& operator+=(const MatrixInterface<MM, ScalarType>& other)
+ {
+ return add_assign(other);
+ }
+
+ template <class MM>
+ derived_type& operator-=(const MatrixInterface<MM, ScalarType>& other)
+ {
+ return subtract_assign(other);
+ }
+
/**
* \brief Checks entries for inf or nan.
* \return false if any entry is inf or nan, else true
@@ -148,6 +178,15 @@ public:
return ret;
} // ... sup_norm(...)
+ derived_type transposed() const
+ {
+ derived_type yy(rows(), cols(), 0.);
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (size_t cc = 0; cc < cols(); ++cc)
+ yy.set_entry(rr, cc, get_entry(cc, rr));
+ return yy;
+ }
+
/**
* \brief Returns the number of entries in the sparsity pattern of the matrix.
*
@@ -243,6 +282,78 @@ public:
} // ... almost_equal(...)
/// \}
+protected:
+ template <class MM>
+ derived_type multiply(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ if (other.rows() != cols())
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, "Dimensions of matrices to be multiplied do not match!");
+ derived_type yy(rows(), other.cols(), 0.);
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (size_t cc = 0; cc < other.cols(); ++cc)
+ for (size_t kk = 0; kk < cols(); ++kk)
+ yy.add_to_entry(rr, cc, get_entry(rr, kk) * other.get_entry(kk, cc));
+ return yy;
+ }
+
+ template <class MM>
+ derived_type add(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ if (other.rows() != rows() || other.cols() != cols())
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, "Dimensions of matrices to be added do not match!");
+ auto new_pattern = pattern() + other.pattern();
+ derived_type yy(rows(), other.cols(), new_pattern);
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (const auto& cc : new_pattern[rr])
+ yy.set_entry(rr, cc, get_entry(rr, cc) + other.get_entry(rr, cc));
+ return yy;
+ }
+
+ template <class MM>
+ derived_type subtract(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ if (other.rows() != rows() || other.cols() != cols())
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, "Dimensions of matrices to be added do not match!");
+ auto new_pattern = pattern() + other.pattern();
+ derived_type yy(rows(), other.cols(), new_pattern);
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (const auto& cc : new_pattern[rr])
+ yy.set_entry(rr, cc, get_entry(rr, cc) - other.get_entry(rr, cc));
+ return yy;
+ }
+
+ template <class MM>
+ derived_type& add_assign(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ if (other.rows() != rows() || other.cols() != cols())
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, "Dimensions of matrices to be added do not match!");
+ const auto this_pattern = pattern();
+ auto new_pattern = this_pattern + other.pattern();
+ if (new_pattern != this_pattern)
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+ "The matrix to be added contains entries that are not in this' pattern!");
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (const auto& cc : this_pattern[rr])
+ add_to_entry(rr, cc, other.get_entry(rr, cc));
+ return this->as_imp();
+ }
+
+ template <class MM>
+ derived_type& subtract_assign(const MatrixInterface<MM, ScalarType>& other) const
+ {
+ if (other.rows() != rows() || other.cols() != cols())
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match, "Dimensions of matrices to be added do not match!");
+ const auto this_pattern = pattern();
+ auto new_pattern = this_pattern + other.pattern();
+ if (new_pattern != this_pattern)
+ DUNE_THROW(XT::Common::Exceptions::shapes_do_not_match,
+ "The matrix to be added contains entries that are not in this' pattern!");
+ for (size_t rr = 0; rr < rows(); ++rr)
+ for (const auto& cc : this_pattern[rr])
+ add_to_entry(rr, cc, -other.get_entry(rr, cc));
+ return this->as_imp();
+ }
+
private:
template <class T, class S>
friend std::ostream& operator<<(std::ostream& /*out*/, const MatrixInterface<T, S>& /*matrix*/);
diff --git a/dune/xt/la/container/pattern.cc b/dune/xt/la/container/pattern.cc
index edeb185c93f2931eeac87f2ce84f853e935bed96..a9ea0030169d9bcde457317a4495991af652d6d0 100644
--- a/dune/xt/la/container/pattern.cc
+++ b/dune/xt/la/container/pattern.cc
@@ -63,6 +63,17 @@ bool SparsityPatternDefault::operator!=(const SparsityPatternDefault& other) con
return vector_of_vectors_ != other.vector_of_vectors_;
}
+SparsityPatternDefault operator+(const SparsityPatternDefault& other) const
+{
+ assert(other.size() == size() && "SparsityPatterns must have the same number of rows for addition!");
+ SparsityPatternDefault ret = *this;
+ for (size_t rr = 0; rr < size(); ++rr)
+ for (const auto& cc : vector_of_vectors_[rr])
+ ret.insert(rr, cc);
+ ret.sort();
+ return ret;
+}
+
void SparsityPatternDefault::insert(const size_t outer_index, const size_t inner_index)
{
assert(outer_index < size() && "Wrong index requested!");
diff --git a/dune/xt/la/container/pattern.hh b/dune/xt/la/container/pattern.hh
index 0a0e84873595ed90358f5ed17ad05eb8788a9d8f..aa36343a3eaaae6988c503d522758e028c420b75 100644
--- a/dune/xt/la/container/pattern.hh
+++ b/dune/xt/la/container/pattern.hh
@@ -45,6 +45,8 @@ public:
bool operator!=(const SparsityPatternDefault& other) const;
+ SparsityPatternDefault operator+(const SparsityPatternDefault& other) const;
+
void insert(const size_t outer_index, const size_t inner_index);
void sort(const size_t outer_index);