From e9cd131056e56d5b3157eaa51d28bcda9f93a9a3 Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Tue, 29 Jan 2019 13:31:19 +0100
Subject: [PATCH] [test.container_matrix] actually test all matrices, add tests
 for clear_row/col, unit_row/col

---
 dune/xt/la/container/common/matrix/dense.hh |  2 +-
 dune/xt/la/container/eigen/dense.hh         | 32 +++++++++
 dune/xt/la/container/eigen/sparse.hh        | 28 +++++++-
 dune/xt/la/container/istl.hh                | 28 +++++++-
 dune/xt/la/container/matrix-interface.hh    |  6 ++
 dune/xt/la/test/container_matrix.py         |  4 +-
 dune/xt/la/test/container_matrix.tpl        | 76 +++++++++++++++++++++
 dune/xt/la/test/matrices.py                 |  1 -
 8 files changed, 171 insertions(+), 6 deletions(-)

diff --git a/dune/xt/la/container/common/matrix/dense.hh b/dune/xt/la/container/common/matrix/dense.hh
index aa6f45f1e..376924c07 100644
--- a/dune/xt/la/container/common/matrix/dense.hh
+++ b/dune/xt/la/container/common/matrix/dense.hh
@@ -381,7 +381,7 @@ public:
         tmp_vec[cc] += vec_entries[ii] * backend_->get_entry_ref(rr, cc);
     }
     for (size_t cc = 0; cc < cols(); ++cc)
-      if (XT::Common::FloatCmp::ne(tmp_vec[cc], 0.))
+      if (XT::Common::FloatCmp::ne(tmp_vec[cc], ScalarType(0.)))
         yy.set_new_entry(cc, tmp_vec[cc]);
   } // void mtv(...)
 
diff --git a/dune/xt/la/container/eigen/dense.hh b/dune/xt/la/container/eigen/dense.hh
index cfac2009e..080562fff 100644
--- a/dune/xt/la/container/eigen/dense.hh
+++ b/dune/xt/la/container/eigen/dense.hh
@@ -565,6 +565,38 @@ public:
     yy.backend().transpose() = backend() * xx.backend();
   }
 
+  template <class T1, class T2>
+  inline void mv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    EigenDenseVector<ScalarType> xx_eigen(xx.size()), yy_eigen(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_eigen[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_eigen[ii] = yy[ii];
+    mv(xx_eigen, yy_eigen);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_eigen[ii]);
+  }
+
+  template <class T1, class T2>
+  inline void mtv(const EigenBaseVector<T1, ScalarType>& xx, EigenBaseVector<T2, ScalarType>& yy) const
+  {
+    yy.backend().transpose() = backend().transpose() * xx.backend();
+  }
+
+  template <class T1, class T2>
+  inline void mtv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    EigenDenseVector<ScalarType> xx_eigen(xx.size()), yy_eigen(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_eigen[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_eigen[ii] = yy[ii];
+    mtv(xx_eigen, yy_eigen);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_eigen[ii]);
+  }
+
   void add_to_entry(const size_t ii, const size_t jj, const ScalarType& value)
   {
     assert(ii < rows());
diff --git a/dune/xt/la/container/eigen/sparse.hh b/dune/xt/la/container/eigen/sparse.hh
index b5ef3a1b9..25b324037 100644
--- a/dune/xt/la/container/eigen/sparse.hh
+++ b/dune/xt/la/container/eigen/sparse.hh
@@ -208,7 +208,7 @@ public:
   ThisType& operator=(const ThisType& other)
   {
     if (this != &other) {
-      backend_ = *other.backend_;
+      *backend_ = *other.backend_;
       mutexes_ = std::make_unique<MutexesType>(other.mutexes_->size());
     }
     return *this;
@@ -286,12 +286,38 @@ public:
     yy.backend().transpose() = backend() * xx.backend();
   }
 
+  template <class T1, class T2>
+  inline void mv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    EigenDenseVector<ScalarType> xx_eigen(xx.size()), yy_eigen(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_eigen[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_eigen[ii] = yy[ii];
+    mv(xx_eigen, yy_eigen);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_eigen[ii]);
+  }
+
   template <class T1, class T2>
   inline void mtv(const EigenBaseVector<T1, ScalarType>& xx, EigenBaseVector<T2, ScalarType>& yy) const
   {
     yy.backend().transpose() = backend().transpose() * xx.backend();
   }
 
+  template <class T1, class T2>
+  inline void mtv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    EigenDenseVector<ScalarType> xx_eigen(xx.size()), yy_eigen(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_eigen[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_eigen[ii] = yy[ii];
+    mtv(xx_eigen, yy_eigen);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_eigen[ii]);
+  }
+
   void add_to_entry(const size_t ii, const size_t jj, const ScalarType& value)
   {
     assert(these_are_valid_indices(ii, jj));
diff --git a/dune/xt/la/container/istl.hh b/dune/xt/la/container/istl.hh
index 0f29903d7..0d0fd7b59 100644
--- a/dune/xt/la/container/istl.hh
+++ b/dune/xt/la/container/istl.hh
@@ -472,7 +472,7 @@ public:
   ThisType& operator=(const ThisType& other)
   {
     if (this != &other) {
-      backend_ = *other.backend_;
+      *backend_ = *other.backend_;
       mutexes_ = std::make_unique<MutexesType>(other.mutexes_->size());
     }
     return *this;
@@ -549,12 +549,38 @@ public:
     backend().mv(xx.backend(), yy.backend());
   }
 
+  template <class T1, class T2>
+  inline void mv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    IstlDenseVector<ScalarType> xx_istl(xx.size()), yy_istl(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_istl[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_istl[ii] = yy[ii];
+    mv(xx_istl, yy_istl);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_istl[ii]);
+  }
+
   inline void mtv(const IstlDenseVector<ScalarType>& xx, IstlDenseVector<ScalarType>& yy) const
   {
     auto& backend_ref = backend();
     backend_ref.mtv(xx.backend(), yy.backend());
   }
 
+  template <class T1, class T2>
+  inline void mtv(const VectorInterface<T1, ScalarType>& xx, VectorInterface<T2, ScalarType>& yy) const
+  {
+    IstlDenseVector<ScalarType> xx_istl(xx.size()), yy_istl(yy.size());
+    for (size_t ii = 0; ii < xx.size(); ++ii)
+      xx_istl[ii] = xx[ii];
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy_istl[ii] = yy[ii];
+    mtv(xx_istl, yy_istl);
+    for (size_t ii = 0; ii < yy.size(); ++ii)
+      yy.set_entry(ii, yy_istl[ii]);
+  }
+
   void add_to_entry(const size_t ii, const size_t jj, const ScalarType& value)
   {
     assert(these_are_valid_indices(ii, jj));
diff --git a/dune/xt/la/container/matrix-interface.hh b/dune/xt/la/container/matrix-interface.hh
index 0abd11611..7de50637a 100644
--- a/dune/xt/la/container/matrix-interface.hh
+++ b/dune/xt/la/container/matrix-interface.hh
@@ -102,6 +102,12 @@ public:
     CHECK_AND_CALL_CRTP(this->as_imp().mv(xx, yy));
   }
 
+  template <class XX, class YY>
+  inline void mtv(const XX& xx, YY& yy) const
+  {
+    CHECK_AND_CALL_CRTP(this->as_imp().mtv(xx, yy));
+  }
+
   inline void add_to_entry(const size_t ii, const size_t jj, const ScalarType& value)
   {
     CHECK_AND_CALL_CRTP(this->as_imp().add_to_entry(ii, jj, value));
diff --git a/dune/xt/la/test/container_matrix.py b/dune/xt/la/test/container_matrix.py
index b52ac7989..670520fa9 100644
--- a/dune/xt/la/test/container_matrix.py
+++ b/dune/xt/la/test/container_matrix.py
@@ -16,5 +16,5 @@ from matrices import matrices, latype, vectors, fieldtypes, vector_filter
 from dune.xt.codegen import typeid_to_typedef_name as safe_name
 
 testtypes = [(safe_name('{}_{}_{}'.format(*mv,f)), latype(mv[0],f), latype(mv[1],f))
-             for mv,f in product(zip(matrices(cache), vectors(cache)), fieldtypes(cache))
-             if vector_filter(mv[1], f)]
\ No newline at end of file
+             for mv,f in product(product(matrices(cache), vectors(cache)), fieldtypes(cache))
+             if vector_filter(mv[1], f)]
diff --git a/dune/xt/la/test/container_matrix.tpl b/dune/xt/la/test/container_matrix.tpl
index f32dbf84e..30ef4b6ef 100644
--- a/dune/xt/la/test/container_matrix.tpl
+++ b/dune/xt/la/test/container_matrix.tpl
@@ -130,6 +130,16 @@ struct MatrixTest_{{T_NAME}} : public ::testing::Test
     sparse_pattern.inner(1).push_back(1); //|-, -, -, x|
     sparse_pattern.inner(2).push_back(3); //|-, -, -, -|
 
+    PatternType sparse_pattern2(dim);
+    for (size_t jj = 0; jj < dim; ++jj)
+      sparse_pattern2.inner(0).push_back(jj); //|x, x, x, x|
+    sparse_pattern2.inner(1).push_back(0);    //|x, x, -, -|
+    sparse_pattern2.inner(1).push_back(1);    //|x, -, x, -|
+    sparse_pattern2.inner(2).push_back(0);    //|x, -, x, -|
+    sparse_pattern2.inner(2).push_back(2);
+    sparse_pattern2.inner(3).push_back(0);
+    sparse_pattern2.inner(3).push_back(2);
+
     // create test matrizes
     MatrixImp matrix_zeros_dense(dim, dim, dense_pattern); // |0, 0, 0, 0|
     for (size_t ii = 0; ii < dim; ++ii) { // |0, 0, 0, 0|
@@ -163,6 +173,11 @@ struct MatrixTest_{{T_NAME}} : public ::testing::Test
     testmatrix_sparse.set_entry(1, 1, ScalarType(1.5)); //|-,   -,   -,    -|
     testmatrix_sparse.set_entry(2, 3, ScalarType(-0.5));
 
+    MatrixImp testmatrix_sparse2(dim, dim, sparse_pattern2);
+    for (size_t ii = 0; ii < dim; ++ii)
+      for (auto&& jj : sparse_pattern2.inner(ii))
+        testmatrix_sparse2.set_entry(ii, jj, 2.);
+
     // create test vectors
     VectorImp vector_zeros(dim); // [0, 0, 0, 0]
     VectorImp vector_ones(dim, ScalarType(1)); // [1, 1, 1, 1]
@@ -330,6 +345,67 @@ struct MatrixTest_{{T_NAME}} : public ::testing::Test
     for (size_t ii = 0; ii < rows; ++ii) {
         DXTC_EXPECT_FLOAT_EQ(res_mtv[ii], res_mv[ii]);
     }
+    // test clear_row/col, unit_row/col
+    auto testm1 = matrix_ones;
+    testm1 *= 2.;
+    testm1.clear_row(0);
+    testm1.clear_row(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (size_t jj = 0; jj < dim; ++jj)
+        DXTC_EXPECT_FLOAT_EQ(testm1.get_entry(ii, jj), ii == 0 || ii == 2 ? ScalarType(0.) : ScalarType(2.));
+    }
+    testm1 = matrix_ones;
+    testm1 *= 2.;
+    testm1.clear_col(0);
+    testm1.clear_col(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (size_t jj = 0; jj < dim; ++jj)
+        DXTC_EXPECT_FLOAT_EQ(testm1.get_entry(ii, jj), jj == 0 || jj == 2 ? ScalarType(0.) : ScalarType(2.));
+    }
+    testm1 = matrix_ones;
+    testm1 *= 2.;
+    testm1.unit_row(0);
+    testm1.unit_row(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (size_t jj = 0; jj < dim; ++jj)
+        DXTC_EXPECT_FLOAT_EQ(testm1.get_entry(ii, jj), ii == 0 || ii == 2 ? ScalarType(ii == jj) : 2.);
+    }
+    testm1 = matrix_ones;
+    testm1 *= 2.;
+    testm1.unit_col(0);
+    testm1.unit_col(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (size_t jj = 0; jj < dim; ++jj)
+        DXTC_EXPECT_FLOAT_EQ(testm1.get_entry(ii, jj), jj == 0 || jj == 2 ? ScalarType(ii == jj) : 2.);
+    }
+    auto testm2 = testmatrix_sparse2;
+    testm2.clear_row(0);
+    testm2.clear_row(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (auto&& jj : sparse_pattern2.inner(ii))
+        DXTC_EXPECT_FLOAT_EQ(testm2.get_entry(ii, jj), ii == 0 || ii == 2 ? ScalarType(0.) : ScalarType(2.));
+    }
+    testm2 = testmatrix_sparse2;
+    testm2.clear_col(0);
+    testm2.clear_col(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (auto&& jj : sparse_pattern2.inner(ii))
+        DXTC_EXPECT_FLOAT_EQ(testm2.get_entry(ii, jj), jj == 0 || jj == 2 ? ScalarType(0.) : ScalarType(2.));
+    }
+    testm2 = testmatrix_sparse2;
+    testm2.unit_row(0);
+    testm2.unit_row(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (auto&& jj : sparse_pattern2.inner(ii))
+        DXTC_EXPECT_FLOAT_EQ(testm2.get_entry(ii, jj), ii == 0 || ii == 2 ? ScalarType(ii == jj) : 2.);
+    }
+    testm2 = testmatrix_sparse2;
+    testm2.unit_col(0);
+    testm2.unit_col(2);
+    for (size_t ii = 0; ii < dim; ++ii) {
+      for (auto&& jj : sparse_pattern2.inner(ii))
+        DXTC_EXPECT_FLOAT_EQ(testm2.get_entry(ii, jj), jj == 0 || jj == 2 ? ScalarType(ii == jj) : 2.);
+    }
   } // void produces_correct_results() const
 }; // struct MatrixTest
 
diff --git a/dune/xt/la/test/matrices.py b/dune/xt/la/test/matrices.py
index 82c645b1e..69589d0c2 100644
--- a/dune/xt/la/test/matrices.py
+++ b/dune/xt/la/test/matrices.py
@@ -16,7 +16,6 @@ from dune.xt import codegen
 def matrices(cache):
   mat = ['CommonDenseMatrix', 'CommonSparseMatrixCsr', 'CommonSparseMatrixCsc',
     'CommonSparseOrDenseMatrixCsr', 'CommonSparseOrDenseMatrixCsc']
-
   if codegen.have_eigen(cache):
       mat.extend(['EigenRowMajorSparseMatrix', 'EigenDenseMatrix'])
   if codegen.have_istl(cache):
-- 
GitLab