From: Julius Witte <julius.witte@iwr.uni-heidelberg.de>
Date: Thu, 21 Sep 2017 18:44:35 +0000 (+0200)
Subject: Added namespace for spectral_assembly(). Base class vmult() and apply_inverse() opera... 
X-Git-Tag: v9.0.0-rc1~1021^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9103de26bd1302189cbdf4fea1083cea3a98d933;p=dealii.git

Added namespace for spectral_assembly(). Base class vmult() and apply_inverse() operate on ArrayView. Deletion of derived class vmult() and apply_inverse().
---

diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h
index 5ac355f53d..03e0dda50a 100644
--- a/include/deal.II/lac/tensor_product_matrix.h
+++ b/include/deal.II/lac/tensor_product_matrix.h
@@ -17,6 +17,7 @@
 #define dealii_tensor_product_matrix_h
 
 
+#include <deal.II/base/array_view.h>
 #include <deal.II/base/config.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/lac/lapack_full_matrix.h>
@@ -83,20 +84,20 @@ public:
   /**
    * Implements a matrix-vector product with the underlying matrix as
    * described in the main documentation of TensorProductMatrixSymmetricSum.
-   * This function is operating on plain pointers, i.e. no check of
-   * array bounds is possible.
+   * This function is operating on ArrayView to allow checks of
+   * array bounds with respect to @p dst and @p src.
    */
-  void vmult (Number *dst,
-              const Number *src) const;
+  void vmult (const ArrayView<Number> &dst,
+              const ArrayView<const Number> &src) const;
 
   /**
    * Implements a matrix-vector product with the underlying matrix as
    * described in the main documentation of TensorProductMatrixSymmetricSum.
-   * This function is operating on plain pointers, i.e. no check of
-   * array bounds is possible.
+   * This function is operating on ArrayView to allow checks of
+   * array bounds with respect to @p dst and @p src.
    */
-  void apply_inverse (Number *dst,
-                      const Number *src) const;
+  void apply_inverse (const ArrayView<Number> &dst,
+                      const ArrayView<const Number> &src) const;
 
 protected:
   /**
@@ -273,34 +274,6 @@ public:
   void reinit (const Table<2,Number> &mass_matrix,
                const Table<2,Number> &derivative_matrix);
 
-  /**
-   * Import functions from base class.
-   */
-  using TensorProductMatrixSymmetricSumBase<dim,Number,size>::vmult;
-
-  /**
-   * Import functions from base class.
-   */
-  using TensorProductMatrixSymmetricSumBase<dim,Number,size>::apply_inverse;
-
-  /**
-   * Implements a matrix-vector product with the underlying matrix as
-   * described in the main documentation of this class. Same as
-   * TensorProductMatrixSymmetricSumBase::vmult() but additionally
-   * providing bound checks of @p dst and @p src.
-   */
-  void vmult (Vector<Number> &dst,
-              const Vector<Number> &src) const;
-
-  /**
-   * Implements a matrix-vector product with the underlying matrix as
-   * described in the main documentation of this class. Same as
-   * TensorProductMatrixSymmetricSumBase::apply_inverse() but additionally
-   * providing bound checks of @p dst and @p src.
-   */
-  void apply_inverse (Vector<Number> &dst,
-                      const Vector<Number> &src) const;
-
 private:
   /**
    * A generic implementation of all reinit() functions based on
@@ -320,6 +293,8 @@ private:
 /** @copydoc TensorProductMatrixSymmetricSum<Number>
  * This is the template specialization for @p Number being
  * VectorizedArray<Number>.
+ *
+ * @author Martin Kronbichler and Julius Witte, 2017
  */
 template <int dim, typename Number, int size>
 class TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
@@ -369,34 +344,6 @@ public:
   void reinit (const Table<2,VectorizedArray<Number> > &mass_matrix,
                const Table<2,VectorizedArray<Number> > &derivative_matrix);
 
-  /**
-   * Import functions from base class.
-   */
-  using TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::vmult;
-
-  /**
-   * Import functions from base class.
-   */
-  using TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::apply_inverse;
-
-  /**
-   * Implements a matrix-vector product with the underlying matrix as
-   * described in the main documentation of this class. Same as
-   * TensorProductMatrixSymmetricSumBase::vmult() but additionally
-   * providing bound checks of @p dst and @p src.
-   */
-  void vmult (AlignedVector<VectorizedArray<Number> > &dst,
-              const AlignedVector<VectorizedArray<Number> > &src) const;
-
-  /**
-   * Implements a matrix-vector product with the underlying matrix as
-   * described in the main documentation of this class. Same as
-   * TensorProductMatrixSymmetricSumBase::apply_inverse() but additionally
-   * providing bound checks of @p dst and @p src.
-   */
-  void apply_inverse (AlignedVector<VectorizedArray<Number> > &dst,
-                      const AlignedVector<VectorizedArray<Number> > &src) const;
-
 private:
   /**
    * A generic implementation of all reinit() functions based on
@@ -416,55 +363,57 @@ private:
 
 #ifndef DOXYGEN
 
-namespace
+namespace internal
 {
-  /**
-   * Compute generalized eigenvalues and eigenvectors of the real
-   * generalized symmetric eigenproblem $A v = \lambda M v$. Since we are
-   * operating on plain pointers we require the size of the matrices beforehand.
-   * Note that the data arrays for the eigenvalues and eigenvectors
-   * have to be initialized to a proper size, too. (no check of array bounds
-   * possible)
-   */
-  template <typename Number>
-  void
-  spectral_assembly (const Number *mass_matrix,
-                     const Number *derivative_matrix,
-                     const unsigned int n_rows,
-                     const unsigned int n_cols,
-                     Number *eigenvalues,
-                     Number *eigenvectors)
+  namespace TensorProductMatrix
   {
-    Assert (n_rows == n_cols, ExcNotImplemented());
-
-    auto &&transpose_fill_nm
-      = [](Number *out, const Number *in, const unsigned int n, const unsigned int m)
+    /**
+     * Compute generalized eigenvalues and eigenvectors of the real
+     * generalized symmetric eigenproblem $A v = \lambda M v$. Since we are
+     * operating on plain pointers we require the size of the matrices beforehand.
+     * Note that the data arrays for the eigenvalues and eigenvectors
+     * have to be initialized to a proper size, too. (no check of array bounds
+     * possible)
+     */
+    template <typename Number>
+    void
+    spectral_assembly (const Number *mass_matrix,
+                       const Number *derivative_matrix,
+                       const unsigned int n_rows,
+                       const unsigned int n_cols,
+                       Number *eigenvalues,
+                       Number *eigenvectors)
     {
-      for (unsigned int mm = 0; mm < m; ++mm)
-        for (unsigned int nn = 0; nn < n; ++nn)
-          out[mm+nn*m] = *(in++);
-    };
-
-    std::vector<Vector<Number> > eigenvecs(n_rows);
-    LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
-    LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
-
-    transpose_fill_nm (&(mass_copy(0,0)), mass_matrix, n_rows, n_cols);
-    transpose_fill_nm (&(deriv_copy(0,0)), derivative_matrix, n_rows, n_cols);
-
-    deriv_copy.compute_generalized_eigenvalues_symmetric (mass_copy, eigenvecs);
-    AssertDimension (eigenvecs.size(), n_rows);
-    for (unsigned int i=0; i<n_rows; ++i)
-      for (unsigned int j=0; j<n_cols; ++j, ++eigenvectors)
-        *eigenvectors = eigenvecs[j][i];
-
-    for (unsigned int i=0; i<n_rows; ++i, ++eigenvalues)
-      *eigenvalues = deriv_copy.eigenvalue(i).real();
+      Assert (n_rows == n_cols, ExcNotImplemented());
+
+      auto &&transpose_fill_nm
+        = [](Number *out, const Number *in, const unsigned int n, const unsigned int m)
+      {
+        for (unsigned int mm = 0; mm < m; ++mm)
+          for (unsigned int nn = 0; nn < n; ++nn)
+            out[mm+nn*m] = *(in++);
+      };
+
+      std::vector<dealii::Vector<Number> > eigenvecs(n_rows);
+      LAPACKFullMatrix<Number> mass_copy(n_rows, n_cols);
+      LAPACKFullMatrix<Number> deriv_copy(n_rows, n_cols);
+
+      transpose_fill_nm (&(mass_copy(0,0)), mass_matrix, n_rows, n_cols);
+      transpose_fill_nm (&(deriv_copy(0,0)), derivative_matrix, n_rows, n_cols);
+
+      deriv_copy.compute_generalized_eigenvalues_symmetric (mass_copy, eigenvecs);
+      AssertDimension (eigenvecs.size(), n_rows);
+      for (unsigned int i=0; i<n_rows; ++i)
+        for (unsigned int j=0; j<n_cols; ++j, ++eigenvectors)
+          *eigenvectors = eigenvecs[j][i];
+
+      for (unsigned int i=0; i<n_rows; ++i, ++eigenvalues)
+        *eigenvalues = deriv_copy.eigenvalue(i).real();
+    }
   }
 }
 
 
-
 template <int dim, typename Number, int size>
 inline
 unsigned int
@@ -495,9 +444,11 @@ template <int dim, typename Number, int size>
 inline
 void
 TensorProductMatrixSymmetricSumBase<dim,Number,size>
-::vmult(Number *dst,
-        const Number *src) const
+::vmult (const ArrayView<Number> &dst_view,
+         const ArrayView<const Number> &src_view) const
 {
+  AssertDimension (dst_view.size(), this->m());
+  AssertDimension (src_view.size(), this->n());
   Threads::Mutex::ScopedLock lock(this->mutex);
   const unsigned int n = Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues[0].size());
   tmp_array.resize_fast(n*2);
@@ -506,6 +457,8 @@ TensorProductMatrixSymmetricSumBase<dim,Number,size>
   eval(AlignedVector<Number> {}, AlignedVector<Number> {},
        AlignedVector<Number> {}, mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows());
   Number *t = tmp_array.begin();
+  const Number *src = src_view.begin();
+  Number *dst = &(dst_view[0]);
 
   if (dim == 1)
     {
@@ -552,9 +505,11 @@ template <int dim, typename Number, int size>
 inline
 void
 TensorProductMatrixSymmetricSumBase<dim,Number,size>
-::apply_inverse(Number *dst,
-                const Number *src) const
+::apply_inverse (const ArrayView<Number> &dst_view,
+                 const ArrayView<const Number> &src_view) const
 {
+  AssertDimension (dst_view.size(), this->n());
+  AssertDimension (src_view.size(), this->m());
   Threads::Mutex::ScopedLock lock(this->mutex);
   const unsigned int n = size > 0 ? size : eigenvalues[0].size();
   tmp_array.resize_fast (Utilities::fixed_power<dim>(n));
@@ -563,6 +518,8 @@ TensorProductMatrixSymmetricSumBase<dim,Number,size>
   eval(AlignedVector<Number>(), AlignedVector<Number>(),
        AlignedVector<Number>(), mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows());
   Number *t = tmp_array.begin();
+  const Number *src = src_view.begin();
+  Number *dst = &(dst_view[0]);
 
   // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
   //       --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
@@ -671,12 +628,13 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
 
       this->eigenvectors[dir].reinit (mass_matrices[dir].n_cols(), mass_matrices[dir].n_rows());
       this->eigenvalues[dir].resize (mass_matrices[dir].n_cols());
-      spectral_assembly<Number> (&(mass_matrices[dir](0,0)),
-                                 &(derivative_matrices[dir](0,0)),
-                                 mass_matrices[dir].n_rows(),
-                                 mass_matrices[dir].n_cols(),
-                                 this->eigenvalues[dir].begin(),
-                                 &(this->eigenvectors[dir](0,0)));
+      internal::TensorProductMatrix
+      ::spectral_assembly<Number> (&(mass_matrices[dir](0,0)),
+                                   &(derivative_matrices[dir](0,0)),
+                                   mass_matrices[dir].n_rows(),
+                                   mass_matrices[dir].n_cols(),
+                                   this->eigenvalues[dir].begin(),
+                                   &(this->eigenvectors[dir](0,0)));
     }
 }
 
@@ -732,34 +690,6 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
 
 
 
-template <int dim, typename Number, int size>
-inline
-void
-TensorProductMatrixSymmetricSum<dim,Number,size>
-::vmult (Vector<Number> &dst,
-         const Vector<Number> &src) const
-{
-  AssertDimension(dst.size(), this->m());
-  AssertDimension(src.size(), this->n());
-  TensorProductMatrixSymmetricSumBase<dim,Number,size>::vmult (dst.begin(), src.begin());
-}
-
-
-
-template <int dim, typename Number, int size>
-inline
-void
-TensorProductMatrixSymmetricSum<dim,Number,size>
-::apply_inverse (Vector<Number> &dst,
-                 const Vector<Number> &src) const
-{
-  AssertDimension (dst.size(), this->n());
-  AssertDimension (src.size(), this->m());
-  TensorProductMatrixSymmetricSumBase<dim,Number,size>::apply_inverse (dst.begin(), src.begin());
-}
-
-
-
 // ------------------------------ vectorized spec.: TensorProductMatrixSymmetricSum   ------------------------------
 
 template <int dim, typename Number, int size>
@@ -849,8 +779,9 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
       Number *eigenvec_begin = eigenvectors_flat.data();
       Number *eigenval_begin = eigenvalues_flat.data();
       for (unsigned int lane = 0; lane < macro_size; ++lane)
-        spectral_assembly<Number> (mass_cbegin+nm*lane, deriv_cbegin+nm*lane, n_rows, n_cols,
-                                   eigenval_begin+n_rows*lane, eigenvec_begin+nm*lane);
+        internal::TensorProductMatrix
+        ::spectral_assembly<Number> (mass_cbegin+nm*lane, deriv_cbegin+nm*lane, n_rows, n_cols,
+                                     eigenval_begin+n_rows*lane, eigenvec_begin+nm*lane);
 
       this->eigenvalues[dir].resize (n_rows);
       this->eigenvectors[dir].reinit (n_rows, n_cols);
@@ -895,34 +826,6 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 
 
 
-template <int dim, typename Number, int size>
-inline
-void
-TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
-::vmult (AlignedVector<VectorizedArray<Number> > &dst,
-         const AlignedVector<VectorizedArray<Number> > &src) const
-{
-  AssertDimension(dst.size(), this->m());
-  AssertDimension(src.size(), this->n());
-  TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::vmult (dst.begin(), src.begin());
-}
-
-
-
-template <int dim, typename Number, int size>
-inline
-void
-TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
-::apply_inverse (AlignedVector<VectorizedArray<Number> > &dst,
-                 const AlignedVector<VectorizedArray<Number> > &src) const
-{
-  AssertDimension (dst.size(), this->n());
-  AssertDimension (src.size(), this->m());
-  TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::apply_inverse (dst.begin(), src.begin());
-}
-
-
-
 #endif
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index 653d9d0c65..678348928d 100644
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -1077,7 +1077,9 @@ namespace internal
             else
               Assert(false, ExcNotImplemented());
 
-            tensor_product_matrix.apply_inverse(&rl1[0], &tmp_rhs[0]);
+            const ArrayView<double> rl1_view {&(rl1[0]), n_inner};
+            const ArrayView<const double> tmp_rhs_view {&(tmp_rhs[0]), n_inner};
+            tensor_product_matrix.apply_inverse(rl1_view, tmp_rhs_view);
 
             for (unsigned int i=0; i<n_inner; ++i)
               lvs(i,k) = -rl1[i];
diff --git a/tests/lac/tensor_product_matrix_01.cc b/tests/lac/tensor_product_matrix_01.cc
index cf4f6c50f0..0ba5436f84 100644
--- a/tests/lac/tensor_product_matrix_01.cc
+++ b/tests/lac/tensor_product_matrix_01.cc
@@ -48,8 +48,11 @@ void do_test(const unsigned int size)
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<double> view1(v1.begin(), v1.size()) ;
+  const ArrayView<double> view2(v2.begin(), v2.size()) ;
+  const ArrayView<double> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
 
@@ -73,7 +76,7 @@ void do_test(const unsigned int size)
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
 }
diff --git a/tests/lac/tensor_product_matrix_02.cc b/tests/lac/tensor_product_matrix_02.cc
index bf7a915d08..2e034bcdd2 100644
--- a/tests/lac/tensor_product_matrix_02.cc
+++ b/tests/lac/tensor_product_matrix_02.cc
@@ -48,8 +48,11 @@ void do_test()
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<double> view1(v1.begin(), v1.size()) ;
+  const ArrayView<double> view2(v2.begin(), v2.size()) ;
+  const ArrayView<double> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
 
@@ -73,7 +76,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
 }
diff --git a/tests/lac/tensor_product_matrix_03.cc b/tests/lac/tensor_product_matrix_03.cc
index ae56c2a463..d82f1f101f 100644
--- a/tests/lac/tensor_product_matrix_03.cc
+++ b/tests/lac/tensor_product_matrix_03.cc
@@ -48,8 +48,11 @@ void do_test()
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<float> view1(v1.begin(), v1.size()) ;
+  const ArrayView<float> view2(v2.begin(), v2.size()) ;
+  const ArrayView<float> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   // add tolerance to account for different BLAS/LAPACK combinations. Float is
   // too inaccurate so numdiff does not work...
@@ -79,7 +82,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   norm = v3.linfty_norm();
   deallog << "Verification of inverse: " << (norm < 5e-3 ? 0. : norm) << std::endl;
diff --git a/tests/lac/tensor_product_matrix_04.cc b/tests/lac/tensor_product_matrix_04.cc
index d2292718e0..840d3aa9bc 100644
--- a/tests/lac/tensor_product_matrix_04.cc
+++ b/tests/lac/tensor_product_matrix_04.cc
@@ -63,8 +63,11 @@ void do_test(const unsigned int size)
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<double> view1(v1.begin(), v1.size()) ;
+  const ArrayView<double> view2(v2.begin(), v2.size()) ;
+  const ArrayView<double> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
 
@@ -91,7 +94,7 @@ void do_test(const unsigned int size)
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
 }
diff --git a/tests/lac/tensor_product_matrix_05.cc b/tests/lac/tensor_product_matrix_05.cc
index 11c0a01bda..1476e76685 100644
--- a/tests/lac/tensor_product_matrix_05.cc
+++ b/tests/lac/tensor_product_matrix_05.cc
@@ -63,8 +63,11 @@ void do_test()
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<double> view1(v1.begin(), v1.size()) ;
+  const ArrayView<double> view2(v2.begin(), v2.size()) ;
+  const ArrayView<double> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
 
@@ -91,7 +94,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
 }
diff --git a/tests/lac/tensor_product_matrix_06.cc b/tests/lac/tensor_product_matrix_06.cc
index f164e422e0..aa320f8054 100644
--- a/tests/lac/tensor_product_matrix_06.cc
+++ b/tests/lac/tensor_product_matrix_06.cc
@@ -63,8 +63,11 @@ void do_test()
   for (unsigned int i=0; i<v1.size(); ++i)
     v1(i) = (2*i+1)%23;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<float> view1(v1.begin(), v1.size()) ;
+  const ArrayView<float> view2(v2.begin(), v2.size()) ;
+  const ArrayView<float> view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   v3 -= v1;
   // add tolerance to account for different BLAS/LAPACK combinations. Float is
   // too inaccurate so numdiff does not work...
@@ -97,7 +100,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(v3, v1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   v3 -= v2;
   norm = v3.linfty_norm();
   deallog << "Verification of inverse: " << (norm < 5e-3 ? 0. : norm) << std::endl;
diff --git a/tests/lac/tensor_product_matrix_vectorized_01.cc b/tests/lac/tensor_product_matrix_vectorized_01.cc
index f5de48716f..15b0c0bea0 100644
--- a/tests/lac/tensor_product_matrix_vectorized_01.cc
+++ b/tests/lac/tensor_product_matrix_vectorized_01.cc
@@ -90,8 +90,11 @@ void do_test(const unsigned int size)
     }) ;
   } ;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<VectorizedArray<double> > view1(v1.begin(), v1.size()) ;
+  const ArrayView<VectorizedArray<double> > view2(v2.begin(), v2.size()) ;
+  const ArrayView<VectorizedArray<double> > view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   subtract_and_assign (v3, v1) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
   deallog << "Verification of vmult and inverse: " << vec_flat.linfty_norm() << std::endl;
@@ -121,7 +124,7 @@ void do_test(const unsigned int size)
 
   full.gauss_jordan();
   full.vmult(w2, w1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   convert_to_vectorized (w2, v3) ;
   subtract_and_assign (v3, v2) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
diff --git a/tests/lac/tensor_product_matrix_vectorized_02.cc b/tests/lac/tensor_product_matrix_vectorized_02.cc
index 0d81df431b..5bd47ce125 100644
--- a/tests/lac/tensor_product_matrix_vectorized_02.cc
+++ b/tests/lac/tensor_product_matrix_vectorized_02.cc
@@ -90,8 +90,11 @@ void do_test()
     }) ;
   } ;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<VectorizedArray<double> > view1(v1.begin(), v1.size()) ;
+  const ArrayView<VectorizedArray<double> > view2(v2.begin(), v2.size()) ;
+  const ArrayView<VectorizedArray<double> > view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   subtract_and_assign (v3, v1) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
   deallog << "Verification of vmult and inverse: " << vec_flat.linfty_norm() << std::endl;
@@ -121,7 +124,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(w2, w1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   convert_to_vectorized (w2, v3) ;
   subtract_and_assign (v3, v2) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
diff --git a/tests/lac/tensor_product_matrix_vectorized_03.cc b/tests/lac/tensor_product_matrix_vectorized_03.cc
index f17a136d42..c465129ad4 100644
--- a/tests/lac/tensor_product_matrix_vectorized_03.cc
+++ b/tests/lac/tensor_product_matrix_vectorized_03.cc
@@ -90,11 +90,14 @@ void do_test()
     }) ;
   } ;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<VectorizedArray<float> > view1(v1.begin(), v1.size()) ;
+  const ArrayView<VectorizedArray<float> > view2(v2.begin(), v2.size()) ;
+  const ArrayView<VectorizedArray<float> > view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   subtract_and_assign (v3, v1) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
-  float norm = vec_flat.linfty_norm();
+  double norm = vec_flat.linfty_norm();
   deallog << "Verification of vmult and inverse: "
           << (norm < 1e-3 ? 0. : norm)  << std::endl;
 
@@ -124,7 +127,7 @@ void do_test()
 
   full.gauss_jordan();
   full.vmult(w2, w1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   convert_to_vectorized (w2, v3) ;
   subtract_and_assign (v3, v2) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
diff --git a/tests/lac/tensor_product_matrix_vectorized_04.cc b/tests/lac/tensor_product_matrix_vectorized_04.cc
index 5884baed3b..a747128bb7 100644
--- a/tests/lac/tensor_product_matrix_vectorized_04.cc
+++ b/tests/lac/tensor_product_matrix_vectorized_04.cc
@@ -84,8 +84,11 @@ void do_test(const unsigned int size)
     }) ;
   } ;
 
-  mat.vmult(v2, v1);
-  mat.apply_inverse(v3, v2);
+  const ArrayView<VectorizedArray<double> > view1(v1.begin(), v1.size()) ;
+  const ArrayView<VectorizedArray<double> > view2(v2.begin(), v2.size()) ;
+  const ArrayView<VectorizedArray<double> > view3(v3.begin(), v3.size()) ;
+  mat.vmult(view2, view1);
+  mat.apply_inverse(view3, view2);
   subtract_and_assign (v3, v1) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
   deallog << "Verification of vmult and inverse: " << vec_flat.linfty_norm() << std::endl;
@@ -115,7 +118,7 @@ void do_test(const unsigned int size)
 
   full.gauss_jordan();
   full.vmult(w2, w1);
-  mat.apply_inverse(v2, v1);
+  mat.apply_inverse(view2, view1);
   convert_to_vectorized (w2, v3) ;
   subtract_and_assign (v3, v2) ;
   vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;