]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added tests for all kinds or reinit() and vectorized template specialization of Tenso...
authorJulius Witte <julius.witte@iwr.uni-heidelberg.de>
Wed, 13 Sep 2017 14:28:39 +0000 (16:28 +0200)
committerJulius Witte <julius.witte@iwr.uni-heidelberg.de>
Fri, 22 Sep 2017 11:31:55 +0000 (13:31 +0200)
15 files changed:
include/deal.II/lac/tensor_product_matrix.h
tests/lac/tensor_product_matrix_04.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_04.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_05.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_05.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_06.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_06.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_01.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_01.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_02.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_02.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_03.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_03.with_lapack=true.output [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_04.cc [new file with mode: 0644]
tests/lac/tensor_product_matrix_vectorized_04.with_lapack=true.output [new file with mode: 0644]

index 8f09f9dcb640662ef139a97f07a52abbdcb8cf5e..99af77b03e74699d2922179c8dec3a4225ccbf23 100644 (file)
@@ -120,7 +120,7 @@ public:
    */
   void vmult (Number *dst,
               const Number *src) const;
-  
+
   /**
    * Implements a matrix-vector product with the underlying matrix as
    * described in the main documentation of this class. Same as the other
@@ -129,33 +129,13 @@ public:
    */
   void apply_inverse (Number *dst,
                       const Number *src) const;
-  
+
 protected:
   /**
    * Constructor.
    */
   TensorProductMatrixSymmetricSumBase () = default ;
 
-  /**
-   * Constructor that is equivalent to the previous constructor and
-   * immediately calling reinit().
-   */
-  TensorProductMatrixSymmetricSumBase (const std::array<Table<2,Number>,dim> &mass_matrix,
-                                       const std::array<Table<2,Number>,dim> &derivative_matrix) ;
-
-  /**
-   * Initializes the matrix to the given mass matrix $M$ and derivative matrix
-   * $A$. Note that the current implementation requires $M$ to be symmetric
-   * and positive definite and $A$ to be symmetric and invertible but not
-   * necessarily positive defininte.
-   */
-  template <typename MatrixArray, typename EigenvalueType, typename EigenvectorType>
-  void
-  fill_data (MatrixArray&& mass_matrices,
-            MatrixArray&& derivative_matrices,
-            EigenvalueType&& eigenvalues,
-            EigenvectorType&& eigenvectors) ;
-
   /**
    * A copy of the @p mass_matrix object passed to the reinit() method.
    */
@@ -176,9 +156,10 @@ protected:
    */
   std::array<Table<2,Number>,dim> eigenvectors;
 
+private:
   /**
-   * An array for temporary data.
-   */
+  * An array for temporary data.
+  */
   mutable AlignedVector<Number> tmp_array;
 
   /**
@@ -208,28 +189,28 @@ public:
    * immediately calling the corresponding reinit().
    */
   TensorProductMatrixSymmetricSum (const std::array<Table<2,Number>, dim> &mass_matrix,
-                                  const std::array<Table<2,Number>, dim> &derivative_matrix) ;
-  
+                                   const std::array<Table<2,Number>, dim> &derivative_matrix) ;
+
   /**
    * Constructor that is equivalent to the first constructor and
    * immediately calling the corresponding reinit().
    */
   TensorProductMatrixSymmetricSum (const std::array<FullMatrix<Number>,dim> &mass_matrix,
                                    const std::array<FullMatrix<Number>,dim> &derivative_matrix) ;
-  
+
   /**
    * Constructor that is equivalent to the first constructor and
-   * immediately calling the corresponding reinit(). 
+   * immediately calling the corresponding reinit().
    */
   TensorProductMatrixSymmetricSum (const FullMatrix<Number> &mass_matrix,
-                                  const FullMatrix<Number> &derivative_matrix) ;
+                                   const FullMatrix<Number> &derivative_matrix) ;
 
   /**
    * Initializes the tensor product matrix to the given mass matrices $M_0,\ldots,M_{dim}$
    * and derivative matrices $A_0,\ldots,A_{dim}$.
    * Note that the current implementation requires each $M_{d}$ to be symmetric
    * and positive definite and every $A_{d}$ to be symmetric and invertible but not
-   * necessarily positive defininte. 
+   * necessarily positive defininte.
    */
   void reinit (const std::array<Table<2,Number>,dim> &mass_matrix,
                const std::array<Table<2,Number>,dim> &derivative_matrix) ;
@@ -240,7 +221,7 @@ public:
    */
   void reinit (const std::array<FullMatrix<Number>,dim> &mass_matrix,
                const std::array<FullMatrix<Number>,dim> &derivative_matrix) ;
-  
+
   /**
    * Initializes the same mass matrix $M$ and derivative matrix $A$ to the given array
    * of mass matrices and array of derivative matrices, respectively.
@@ -277,14 +258,14 @@ public:
 
 private:
   /**
-   * A generic implementation of all reinit() functions based on 
+   * A generic implementation of all reinit() functions based on
    * perfect forwarding, that makes it possible to pass lvalue as well
    * as rvalue arguments. MatrixArray has to be convertible to the underlying
    * type of the bass class' members mass_matrices and derivative_matrices.
    */
   template <typename MatrixArray>
   void reinit_impl (MatrixArray &&mass_matrix,
-                   MatrixArray &&derivative_matrix) ;
+                    MatrixArray &&derivative_matrix) ;
 };
 
 
@@ -311,17 +292,17 @@ public:
 
   /**
    * Constructor that is equivalent to the first constructor and
-   * immediately calling the corresponding reinit(). 
+   * immediately calling the corresponding reinit().
    */
   TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray<Number> > &mass_matrix,
-                                  const Table<2,VectorizedArray<Number> > &derivative_matrix) ;
+                                   const Table<2,VectorizedArray<Number> > &derivative_matrix) ;
 
   /**
    * Initializes the tensor product matrix to the given mass matrices $M_0,\ldots,M_{dim}$
    * and derivative matrices $A_0,\ldots,A_{dim}$.
    * Note that the current implementation requires each $M_{d}$ to be symmetric
    * and positive definite and every $A_{d}$ to be symmetric and invertible but not
-   * necessarily positive defininte. 
+   * necessarily positive defininte.
    */
   void reinit (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
                const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix) ;
@@ -352,14 +333,14 @@ public:
 
 private:
   /**
-   * A generic implementation of all reinit() functions based on 
+   * A generic implementation of all reinit() functions based on
    * perfect forwarding, that makes it possible to pass lvalue as well
    * as rvalue arguments. MatrixArray has to be convertible to the underlying
    * type of the bass class' members mass_matrices and derivative_matrices.
    */
   template <typename MatrixArray>
   void reinit_impl (MatrixArray &&mass_matrix,
-                   MatrixArray &&derivative_matrix) ;
+                    MatrixArray &&derivative_matrix) ;
 };
 
 
@@ -378,12 +359,13 @@ namespace
    * 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)
+  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)
   {
     Assert (n_rows == n_cols, ExcNotImplemented()) ;
 
@@ -415,28 +397,6 @@ namespace
 
 
 
-template <int dim, typename Number, int size>
-template <typename MatrixArray, typename EigenvalueType, typename EigenvectorType>
-inline
-void
-TensorProductMatrixSymmetricSumBase<dim,Number,size>
-::fill_data (MatrixArray&& mass_matrices,
-            MatrixArray&& derivative_matrices,
-            EigenvalueType&& eigenvalues,
-            EigenvectorType&& eigenvectors)
-{
-  AssertDimension (mass_matrices.size(), dim) ;
-  AssertDimension (eigenvalues.size(), dim) ;
-  AssertDimension (eigenvectors.size(), dim) ;
-
-  this->mass_matrix = std::forward<MatrixArray>(mass_matrices) ;
-  this->derivative_matrix = std::forward<MatrixArray>(derivative_matrices) ;
-  this->eigenvalues = std::forward<EigenvalueType>(eigenvalues) ;
-  this->eigenvectors = std::forward<EigenvectorType>(eigenvectors) ;
-}
-
-
-
 template <int dim, typename Number, int size>
 inline
 unsigned int
@@ -600,9 +560,9 @@ template <int dim, typename Number, int size>
 inline
 TensorProductMatrixSymmetricSum<dim,Number,size>
 ::TensorProductMatrixSymmetricSum (const std::array<Table<2,Number>, dim> &mass_matrix,
-                                  const std::array<Table<2,Number>, dim> &derivative_matrix)
+                                   const std::array<Table<2,Number>, dim> &derivative_matrix)
 {
-  reinit_impl (mass_matrix, derivative_matrix) ;
+  reinit (mass_matrix, derivative_matrix) ;
 }
 
 
@@ -635,13 +595,13 @@ inline
 void
 TensorProductMatrixSymmetricSum<dim,Number,size>
 ::reinit_impl (MatrixArray &&mass_matrices_,
-              MatrixArray &&derivative_matrices_)
+               MatrixArray &&derivative_matrices_)
 {
   auto &&mass_matrices = std::forward<MatrixArray>(mass_matrices_) ;
   auto &&derivative_matrices = std::forward<MatrixArray>(derivative_matrices_) ;
-  
-  std::array<Table<2,Number>,dim> eigenvectors ;
-  std::array<AlignedVector<Number>, dim> eigenvalues ;
+  this->mass_matrix = mass_matrices ;
+  this->derivative_matrix = derivative_matrices ;
+
   for (int dir = 0; dir < dim; ++dir)
     {
       Assert (size == -1 || (size > 0 && static_cast<unsigned int>(size) == mass_matrices[dir].n_rows()),
@@ -650,19 +610,15 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
       AssertDimension (mass_matrices[dir].n_rows(), derivative_matrices[dir].n_rows());
       AssertDimension (mass_matrices[dir].n_rows(), derivative_matrices[dir].n_cols());
 
-      eigenvectors[dir].reinit (mass_matrices[dir].n_cols(), mass_matrices[dir].n_rows()) ;
-      eigenvalues[dir].resize (mass_matrices[dir].n_cols()) ;
+      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()
-                                 , eigenvalues[dir].begin()
-                                 , &(eigenvectors[dir](0,0))) ;
+                                 , this->eigenvalues[dir].begin()
+                                 , &(this->eigenvectors[dir](0,0))) ;
     }
-
-  TensorProductMatrixSymmetricSumBase<dim,Number,size>
-    ::fill_data (std::forward<MatrixArray>(mass_matrices), std::forward<MatrixArray>(derivative_matrices),
-                std::move(eigenvalues), std::move(eigenvectors)) ;
 }
 
 
@@ -690,9 +646,9 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
   std::array<Table<2,Number>,dim> deriv_copy ;
 
   std::transform (mass_matrix.cbegin(), mass_matrix.cend(), mass_copy.begin(),
-                 [] (const FullMatrix<Number>& m) ->Table<2,Number> {return m;}) ;
+                  [] (const FullMatrix<Number> &m) ->Table<2,Number> {return m;}) ;
   std::transform (derivative_matrix.cbegin(), derivative_matrix.cend(), deriv_copy.begin(),
-                 [] (const FullMatrix<Number>& m) ->Table<2,Number> {return m;}) ;
+                  [] (const FullMatrix<Number> &m) ->Table<2,Number> {return m;}) ;
 
   reinit_impl (std::move(mass_copy), std::move(deriv_copy)) ;
 }
@@ -708,7 +664,7 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
 {
   std::array<Table<2,Number>,dim> mass_matrices ;
   std::array<Table<2,Number>,dim> derivative_matrices ;
-  
+
   std::fill (mass_matrices.begin(), mass_matrices.end(), mass_matrix) ;
   std::fill (derivative_matrices.begin(), derivative_matrices.end(), derivative_matrix) ;
 
@@ -724,8 +680,8 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
 ::vmult (Vector<Number> &dst,
          const Vector<Number> &src) const
 {
-  AssertDimension(dst.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
-  AssertDimension(src.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
+  AssertDimension(dst.size(), this->m()) ;
+  AssertDimension(src.size(), this->n()) ;
   TensorProductMatrixSymmetricSumBase<dim,Number,size>::vmult (dst.begin(), src.begin());
 }
 
@@ -738,63 +694,14 @@ TensorProductMatrixSymmetricSum<dim,Number,size>
 ::apply_inverse (Vector<Number> &dst,
                  const Vector<Number> &src) const
 {
-  AssertDimension (dst.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
-  AssertDimension (src.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
+  AssertDimension (dst.size(), this->n()) ;
+  AssertDimension (src.size(), this->m()) ;
   TensorProductMatrixSymmetricSumBase<dim,Number,size>::apply_inverse (dst.begin(), src.begin());
 }
 
 
 
-// template <int dim, typename Number, int size>
-// inline
-// FullMatrix<Number>
-// TensorProductMatrixSymmetricSum<dim,Number,size>
-// ::get_full_matrix () const
-// {
-//   const auto& mass_matrix = TensorProductMatrixSymmetricSumBase<dim,Number,size>::mass_matrix ;
-//   const auto& derivative_matrix = this->derivative_matrix ;
-//   const auto& eigenvalues = this->eigenvalues ;
-
-//   FullMatrix<Number> matrix {Utilities::fixed_power<dim>(mass_matrix[0].n_rows())} ;
-//   const unsigned int stride = size > 0 ? size : eigenvalues[0].size() ;
-
-//   if (dim == 1)
-//     matrix.Table<2,Number>::fill (&(derivative_matrix[0](0,0)), true) ;
-
-//   else if (dim == 2)
-//     {
-//       for (unsigned int i1 = 0; i1 < stride; ++i1)
-//  for (unsigned int j1 = 0; j1 < stride; ++j1)
-//    for (unsigned int i0 = 0; i0 < stride; ++i0)
-//      for (unsigned int j0 = 0; j0 < stride; ++j0)
-//        matrix(i1*stride+i0, j1*stride+j0)
-//    = mass_matrix[1](i1,j1) * derivative_matrix[0](i0,j0)
-//    + derivative_matrix[1](i1,j1) * mass_matrix[0](i0,j0) ;
-//     }
-
-//   else if (dim == 3)
-//     {
-//       const unsigned int stride2 = stride * stride ;
-//       for (unsigned int i2 = 0; i2 < stride; ++i2)
-//  for (unsigned int j2 = 0; j2 < stride; ++j2)
-//    for (unsigned int i1 = 0; i1 < stride; ++i1)
-//      for (unsigned int j1 = 0; j1 < stride; ++j1)
-//        for (unsigned int i0 = 0; i0 < stride; ++i0)
-//    for (unsigned int j0 = 0; j0 < stride; ++j0)
-//      matrix(i2*stride2+i1*stride+i0, j2*stride2+j1*stride+j0)
-//        = mass_matrix[2](i2,j2) * mass_matrix[1](i1,j1) * derivative_matrix[0](i0,j0)
-//        + mass_matrix[2](i2,j2) * derivative_matrix[1](i1,j1) * mass_matrix[0](i0,j0)
-//        + derivative_matrix[2](i2,j2) * mass_matrix[1](i1,j1) * mass_matrix[0](i0,j0) ;
-//     }
-
-//   else
-//     Assert (false, ExcNotImplemented()) ;
-
-//   return matrix ;
-// }
-
-
-// ------------------------------ vectorized spez.: TensorProductMatrixSymmetricSum   ------------------------------
+// ------------------------------ vectorized spec.: TensorProductMatrixSymmetricSum   ------------------------------
 
 template <int dim, typename Number, int size>
 inline
@@ -811,7 +718,7 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 ::TensorProductMatrixSymmetricSum (const std::array<Table<2,VectorizedArray<Number> >,dim> &mass_matrix,
                                    const std::array<Table<2,VectorizedArray<Number> >,dim> &derivative_matrix)
 {
-  reinit_impl (mass_matrix, derivative_matrix) ;
+  reinit (mass_matrix, derivative_matrix) ;
 }
 
 
@@ -820,15 +727,9 @@ template <int dim, typename Number, int size>
 inline
 TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 ::TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray<Number> > &mass_matrix,
-                                  const Table<2,VectorizedArray<Number> > &derivative_matrix)
+                                   const Table<2,VectorizedArray<Number> > &derivative_matrix)
 {
-  std::array<Table<2,VectorizedArray<Number> >,dim> mass_matrices ;
-  std::array<Table<2,VectorizedArray<Number> >,dim> derivative_matrices ;
-
-  std::fill (mass_matrices.begin(), mass_matrices.end(), mass_matrix) ;
-  std::fill (derivative_matrices.begin(), derivative_matrices.end(), derivative_matrix) ;
-  
-  reinit_impl (mass_matrices, derivative_matrices) ;
+  reinit (mass_matrix, derivative_matrix) ;
 }
 
 
@@ -839,19 +740,33 @@ inline
 void
 TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 ::reinit_impl (MatrixArray &&mass_matrices_,
-              MatrixArray &&derivative_matrices_)
+               MatrixArray &&derivative_matrices_)
 {
   auto &&mass_matrix = std::forward<MatrixArray>(mass_matrices_) ;
   auto &&derivative_matrix = std::forward<MatrixArray>(derivative_matrices_) ;
-  std::array<Table<2,VectorizedArray<Number> >,dim> eigenvectors ;
-  std::array<AlignedVector<VectorizedArray<Number> >, dim> eigenvalues ;
-  
+  this->mass_matrix = mass_matrix ;
+  this->derivative_matrix = derivative_matrix ;
+
   constexpr unsigned int macro_size = VectorizedArray<Number>::n_array_elements ;
+  const unsigned int nm_flat_size
+    = (size > 0)
+      ? (Utilities::fixed_int_power<size,dim>::value
+         * Utilities::fixed_int_power<size,dim>::value * macro_size)
+      : (Utilities::fixed_power<dim>(mass_matrix[0].n_rows())
+         * Utilities::fixed_power<dim>(mass_matrix[0].n_rows()) * macro_size) ;
+  const unsigned int n_flat_size
+    = (size > 0)
+      ? Utilities::fixed_int_power<size,dim>::value * macro_size
+      : Utilities::fixed_power<dim>(mass_matrix[0].n_rows()) * macro_size ;
 
   std::vector<Number> mass_matrix_flat ;
   std::vector<Number> deriv_matrix_flat ;
   std::vector<Number> eigenvalues_flat ;
   std::vector<Number> eigenvectors_flat ;
+  mass_matrix_flat.reserve (nm_flat_size) ;
+  deriv_matrix_flat.reserve (nm_flat_size) ;
+  eigenvalues_flat.reserve (n_flat_size) ;
+  eigenvectors_flat.reserve (nm_flat_size) ;
   std::array<unsigned int,macro_size> offsets_nm ;
   std::array<unsigned int,macro_size> offsets_n ;
   for (int dir = 0; dir < dim; ++dir)
@@ -866,14 +781,15 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
       const unsigned int n_rows = mass_matrix[dir].n_rows() ;
       const unsigned int n_cols = mass_matrix[dir].n_cols() ;
       const unsigned int nm = n_rows * n_cols ;
+
       mass_matrix_flat.resize (macro_size*nm) ;
       deriv_matrix_flat.resize (macro_size*nm) ;
       eigenvalues_flat.resize (macro_size*n_rows) ;
       eigenvectors_flat.resize (macro_size*nm) ;
-      std::generate (offsets_nm.begin(), offsets_nm.end(),
-                     [=, i=unsigned {0}] () mutable {return nm*(i++);}) ;
-      std::generate (offsets_n.begin(), offsets_n.end(),
-                     [=, i=unsigned {0}] () mutable {return n_rows*(i++);}) ;
+      for (unsigned int vv=0; vv<macro_size; ++vv)
+        offsets_nm[vv] = nm * vv ;
+      for (unsigned int vv=0; vv<macro_size; ++vv)
+        offsets_n[vv] = n_rows * vv ;
 
       vectorized_transpose_and_store (false, nm, &(mass_matrix[dir](0,0))
                                       , offsets_nm.cbegin(), mass_matrix_flat.data()) ;
@@ -884,6 +800,7 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
       const Number *deriv_cbegin = deriv_matrix_flat.data() ;
       Number *eigenvec_begin = eigenvectors_flat.data() ;
       Number *eigenval_begin = eigenvalues_flat.data() ;
+
       spectral_assembly<Number> (mass_cbegin, deriv_cbegin, n_rows, n_cols
                                  , eigenval_begin, eigenvec_begin) ;
       for (unsigned int lane = 1; lane < macro_size; ++lane)
@@ -896,17 +813,13 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
                                      , eigenval_begin, eigenvec_begin) ;
         }
 
-      eigenvalues[dir].resize (n_rows) ;
-      eigenvectors[dir].reinit (n_rows, n_cols) ;
+      this->eigenvalues[dir].resize (n_rows) ;
+      this->eigenvectors[dir].reinit (n_rows, n_cols) ;
       vectorized_load_and_transpose (n_rows, eigenvalues_flat.data()
                                      , offsets_n.cbegin(), this->eigenvalues[dir].begin()) ;
       vectorized_load_and_transpose (nm, eigenvectors_flat.data()
                                      , offsets_nm.cbegin(), &(this->eigenvectors[dir](0,0))) ;
     }
-
-  TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>
-    ::fill_data (std::forward<MatrixArray>(mass_matrix), std::forward<MatrixArray>(derivative_matrix),
-                std::move(eigenvalues), std::move(eigenvectors)) ;
 }
 
 
@@ -923,6 +836,24 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 
 
 
+template <int dim, typename Number, int size>
+inline
+void
+TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
+::reinit (const Table<2,VectorizedArray<Number> > &mass_matrix,
+          const Table<2,VectorizedArray<Number> > &derivative_matrix)
+{
+  std::array<Table<2,VectorizedArray<Number> >,dim> mass_matrices ;
+  std::array<Table<2,VectorizedArray<Number> >,dim> derivative_matrices ;
+
+  std::fill (mass_matrices.begin(), mass_matrices.end(), mass_matrix) ;
+  std::fill (derivative_matrices.begin(), derivative_matrices.end(), derivative_matrix) ;
+
+  reinit_impl (std::move(mass_matrices), std::move(derivative_matrices)) ;
+}
+
+
+
 template <int dim, typename Number, int size>
 inline
 void
@@ -930,8 +861,8 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 ::vmult (AlignedVector<VectorizedArray<Number> > &dst,
          const AlignedVector<VectorizedArray<Number> > &src) const
 {
-  AssertDimension(dst.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
-  AssertDimension(src.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
+  AssertDimension(dst.size(), this->m()) ;
+  AssertDimension(src.size(), this->n()) ;
   TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::vmult (dst.begin(), src.begin());
 }
 
@@ -944,8 +875,8 @@ TensorProductMatrixSymmetricSum<dim,VectorizedArray<Number>,size>
 ::apply_inverse (AlignedVector<VectorizedArray<Number> > &dst,
                  const AlignedVector<VectorizedArray<Number> > &src) const
 {
-  AssertDimension (dst.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
-  AssertDimension (src.size(), Utilities::fixed_power<dim>(this->eigenvalues[0].size()));
+  AssertDimension (dst.size(), this->n()) ;
+  AssertDimension (src.size(), this->m()) ;
   TensorProductMatrixSymmetricSumBase<dim,VectorizedArray<Number>,size>::apply_inverse (dst.begin(), src.begin());
 }
 
diff --git a/tests/lac/tensor_product_matrix_04.cc b/tests/lac/tensor_product_matrix_04.cc
new file mode 100644 (file)
index 0000000..4fc4738
--- /dev/null
@@ -0,0 +1,117 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Similar to tensor_product_matrix_01.cc unless testing with
+// different mass and laplace matrices for each tensor direction, respectively.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim>
+void do_test(const unsigned int size)
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<double> init_mass(size, size);
+  FullMatrix<double> init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3. ;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6. ;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6. ;
+      init_laplace(i,i) = 2. ;
+      if (i > 0)
+        init_laplace(i,i-1) = -1. ;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1. ;
+    }
+
+  std::array<FullMatrix<double>, dim> mass ;
+  std::array<FullMatrix<double>, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= 4./3. ;
+          init_laplace(i,i) *= 5./4. ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+
+  TensorProductMatrixSymmetricSum<dim,double> mat;
+  mat.reinit(mass, laplace);
+  Vector<double> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk);
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk);
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii) * mass[1](j,jj)*mass[0](k,kk)
+                      + mass[2](i,ii) * (laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk));
+  full.vmult(v3, v1);
+  v3 -= v2;
+  deallog << "Verifiction of vmult: " << v3.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1>(1);
+  do_test<1>(2);
+  do_test<1>(5);
+  do_test<2>(1);
+  do_test<2>(2);
+  do_test<2>(5);
+  do_test<2>(11);
+  do_test<3>(1);
+  do_test<3>(2);
+  do_test<3>(3);
+  do_test<3>(7);
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_04.with_lapack=true.output b/tests/lac/tensor_product_matrix_04.with_lapack=true.output
new file mode 100644 (file)
index 0000000..bde376c
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_05.cc b/tests/lac/tensor_product_matrix_05.cc
new file mode 100644 (file)
index 0000000..d0f84f6
--- /dev/null
@@ -0,0 +1,117 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Similar to tensor_product_matrix_02.cc unless testing with
+// different mass and laplace matrices for each tensor direction, respectively.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<double> init_mass(size, size);
+  FullMatrix<double> init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3. ;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6. ;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6. ;
+      init_laplace(i,i) = 2. ;
+      if (i > 0)
+        init_laplace(i,i-1) = -1. ;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1. ;
+    }
+
+  std::array<FullMatrix<double>, dim> mass ;
+  std::array<FullMatrix<double>, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= 4./3. ;
+          init_laplace(i,i) *= 5./4. ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+
+  TensorProductMatrixSymmetricSum<dim,double,size> mat;
+  mat.reinit(mass, laplace);
+  Vector<double> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  deallog << "Verification of vmult and inverse: " << v3.linfty_norm() << std::endl;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk);
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk);
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii) * mass[1](j,jj)*mass[0](k,kk)
+                      + mass[2](i,ii) * (laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk));
+  full.vmult(v3, v1);
+  v3 -= v2;
+  deallog << "Verifiction of vmult: " << v3.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  deallog << "Verification of inverse: " << v3.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_05.with_lapack=true.output b/tests/lac/tensor_product_matrix_05.with_lapack=true.output
new file mode 100644 (file)
index 0000000..bde376c
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verifiction of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_06.cc b/tests/lac/tensor_product_matrix_06.cc
new file mode 100644 (file)
index 0000000..2ec0b53
--- /dev/null
@@ -0,0 +1,124 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Similar to tensor_product_matrix_03.cc unless testing with
+// different mass and laplace matrices for each tensor direction, respectively.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  FullMatrix<float> init_mass(size, size);
+  FullMatrix<float> init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3. ;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6. ;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6. ;
+      init_laplace(i,i) = 2. ;
+      if (i > 0)
+        init_laplace(i,i-1) = -1. ;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1. ;
+    }
+
+  std::array<FullMatrix<float>, dim> mass ;
+  std::array<FullMatrix<float>, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= 4./3. ;
+          init_laplace(i,i) *= 5./4. ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+
+  TensorProductMatrixSymmetricSum<dim,float,size> mat;
+  mat.reinit(mass, laplace);
+  Vector<float> v1(mat.m()), v2(mat.m()), v3(mat.m());
+  for (unsigned int i=0; i<v1.size(); ++i)
+    v1(i) = (2*i+1)%23;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  v3 -= v1;
+  // add tolerance to account for different BLAS/LAPACK combinations. Float is
+  // too inaccurate so numdiff does not work...
+  double norm = v3.linfty_norm();
+  deallog << "Verification of vmult and inverse: "
+          << (norm < 1e-3 ? 0. : norm)  << std::endl;
+
+  FullMatrix<float> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk);
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk);
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii) * mass[1](j,jj)*mass[0](k,kk)
+                      + mass[2](i,ii) * (laplace[1](j,jj)*mass[0](k,kk) + mass[1](j,jj)*laplace[0](k,kk));
+  full.vmult(v3, v1);
+  v3 -= v2;
+
+  norm = v3.linfty_norm();
+  deallog << "Verifiction of vmult: " << (norm < 1e-4 ? 0. : norm) << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(v3, v1);
+  mat.apply_inverse(v2, v1);
+  v3 -= v2;
+  norm = v3.linfty_norm();
+  deallog << "Verification of inverse: " << (norm < 5e-3 ? 0. : norm) << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_06.with_lapack=true.output b/tests/lac/tensor_product_matrix_06.with_lapack=true.output
new file mode 100644 (file)
index 0000000..2fc7805
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verifiction of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
diff --git a/tests/lac/tensor_product_matrix_vectorized_01.cc b/tests/lac/tensor_product_matrix_vectorized_01.cc
new file mode 100644 (file)
index 0000000..ec37cb8
--- /dev/null
@@ -0,0 +1,151 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Same as 'tensor_product_matrix_04.cc' unless that we replaced the scalar data
+// type 'double' by the vectorized data type 'VectorizedArray<double>'.
+// Note, all lanes compute the same.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/base/vectorization.h>
+
+template <int dim>
+void do_test(const unsigned int size)
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  Table<2,VectorizedArray<double> > init_mass(size, size);
+  Table<2,VectorizedArray<double> > init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3.;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6.;
+      init_laplace(i,i) = 2.;
+      if (i > 0)
+        init_laplace(i,i-1) = -1.;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1.;
+    }
+  std::array<Table<2,VectorizedArray<double> >, dim> mass ;
+  std::array<Table<2,VectorizedArray<double> >, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= make_vectorized_array<double>(4./3.) ;
+          init_laplace(i,i) *= make_vectorized_array<double>(5./4.) ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+  TensorProductMatrixSymmetricSum<dim,VectorizedArray<double> > mat;
+  mat.reinit(mass, laplace);
+
+  Vector<double> w1(mat.m()), w2(mat.m()) ;
+  for (unsigned int i=0; i<w1.size(); ++i)
+    w1[i] = (2*i+1)%23 ;
+
+  auto convert_to_vectorized =
+    [](const Vector<double> &in
+       , AlignedVector<VectorizedArray<double> > &out)
+  {
+    std::transform (in.begin(), in.end(), out.begin(),
+                    [](const auto &val)
+    {
+      return make_vectorized_array(val);
+    }) ;
+  } ;
+  AlignedVector<VectorizedArray<double> > v1(w1.size()), v2(w1.size()), v3(w1.size());
+  convert_to_vectorized (w1, v1) ;
+
+  constexpr unsigned int macro_size = VectorizedArray<double>::n_array_elements ;
+  Vector<double> vec_flat(v1.size()*macro_size) ;
+  std::array<unsigned int,macro_size> offsets ;
+  for (unsigned int i=0; i<macro_size; ++i)
+    offsets[i] = v1.size() * i ;
+  auto subtract_and_assign =
+    [](AlignedVector<VectorizedArray<double> > &lhs
+       , const AlignedVector<VectorizedArray<double> > &rhs)
+  {
+    std::transform (lhs.begin(), lhs.end(), rhs.begin(), lhs.begin(),
+                    [](const auto lval, const auto rval)
+    {
+      return lval - rval;
+    }) ;
+  } ;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  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;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk)[0];
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0];
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii)[0] * mass[1](j,jj)[0]*mass[0](k,kk)[0]
+                      + mass[2](i,ii)[0] * (laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0]);
+  full.vmult(w2, w1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of vmult: " << vec_flat.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(w2, w1);
+  mat.apply_inverse(v2, v1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of inverse: " << vec_flat.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1>(1);
+  do_test<1>(2);
+  do_test<1>(5);
+  do_test<2>(1);
+  do_test<2>(2);
+  do_test<2>(5);
+  do_test<2>(11);
+  do_test<3>(1);
+  do_test<3>(2);
+  do_test<3>(3);
+  do_test<3>(7);
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_vectorized_01.with_lapack=true.output b/tests/lac/tensor_product_matrix_vectorized_01.with_lapack=true.output
new file mode 100644 (file)
index 0000000..4ef63ad
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_vectorized_02.cc b/tests/lac/tensor_product_matrix_vectorized_02.cc
new file mode 100644 (file)
index 0000000..db59761
--- /dev/null
@@ -0,0 +1,151 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Same as 'tensor_product_matrix_05' unless that we replaced the scalar data
+// type 'double' by the vectorized data type 'VectorizedArray<double>'.
+// Note, all lanes compute the same.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/base/vectorization.h>
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  Table<2,VectorizedArray<double> > init_mass(size, size);
+  Table<2,VectorizedArray<double> > init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3.;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6.;
+      init_laplace(i,i) = 2.;
+      if (i > 0)
+        init_laplace(i,i-1) = -1.;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1.;
+    }
+  std::array<Table<2,VectorizedArray<double> >, dim> mass ;
+  std::array<Table<2,VectorizedArray<double> >, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= make_vectorized_array<double>(4./3.) ;
+          init_laplace(i,i) *= make_vectorized_array<double>(5./4.) ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+  TensorProductMatrixSymmetricSum<dim,VectorizedArray<double> > mat;
+  mat.reinit(mass, laplace);
+
+  Vector<double> w1(mat.m()), w2(mat.m()) ;
+  for (unsigned int i=0; i<w1.size(); ++i)
+    w1[i] = (2*i+1)%23 ;
+
+  auto convert_to_vectorized =
+    [](const Vector<double> &in
+       , AlignedVector<VectorizedArray<double> > &out)
+  {
+    std::transform (in.begin(), in.end(), out.begin(),
+                    [](const auto &val)
+    {
+      return make_vectorized_array(val);
+    }) ;
+  } ;
+  AlignedVector<VectorizedArray<double> > v1(w1.size()), v2(w1.size()), v3(w1.size());
+  convert_to_vectorized (w1, v1) ;
+
+  constexpr unsigned int macro_size = VectorizedArray<double>::n_array_elements ;
+  Vector<double> vec_flat(v1.size()*macro_size) ;
+  std::array<unsigned int,macro_size> offsets ;
+  for (unsigned int i=0; i<macro_size; ++i)
+    offsets[i] = v1.size() * i ;
+  auto subtract_and_assign =
+    [](AlignedVector<VectorizedArray<double> > &lhs
+       , const AlignedVector<VectorizedArray<double> > &rhs)
+  {
+    std::transform (lhs.begin(), lhs.end(), rhs.begin(), lhs.begin(),
+                    [](const auto lval, const auto rval)
+    {
+      return lval - rval;
+    }) ;
+  } ;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  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;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk)[0];
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0];
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii)[0] * mass[1](j,jj)[0]*mass[0](k,kk)[0]
+                      + mass[2](i,ii)[0] * (laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0]);
+  full.vmult(w2, w1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of vmult: " << vec_flat.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(w2, w1);
+  mat.apply_inverse(v2, v1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of inverse: " << vec_flat.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_vectorized_02.with_lapack=true.output b/tests/lac/tensor_product_matrix_vectorized_02.with_lapack=true.output
new file mode 100644 (file)
index 0000000..4ef63ad
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
diff --git a/tests/lac/tensor_product_matrix_vectorized_03.cc b/tests/lac/tensor_product_matrix_vectorized_03.cc
new file mode 100644 (file)
index 0000000..0e73c3a
--- /dev/null
@@ -0,0 +1,155 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Same as 'tensor_product_matrix_06.cc' unless that we replaced the scalar data
+// type 'float' by the vectorized data type 'VectorizedArray<float>'.
+// Note, all lanes compute the same.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/base/vectorization.h>
+
+template <int dim, int size>
+void do_test()
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  Table<2,VectorizedArray<float> > init_mass(size, size);
+  Table<2,VectorizedArray<float> > init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3.;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6.;
+      init_laplace(i,i) = 2.;
+      if (i > 0)
+        init_laplace(i,i-1) = -1.;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1.;
+    }
+  std::array<Table<2,VectorizedArray<float> >, dim> mass ;
+  std::array<Table<2,VectorizedArray<float> >, dim> laplace ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    {
+      for (unsigned int i=0; i<size; ++i)
+        {
+          init_mass(i,i) *= make_vectorized_array<float>(4./3.) ;
+          init_laplace(i,i) *= make_vectorized_array<float>(5./4.) ;
+        }
+      mass[dir] = init_mass ;
+      laplace[dir] = init_laplace ;
+    }
+  TensorProductMatrixSymmetricSum<dim,VectorizedArray<float> > mat;
+  mat.reinit(mass, laplace);
+
+  Vector<float> w1(mat.m()), w2(mat.m()) ;
+  for (unsigned int i=0; i<w1.size(); ++i)
+    w1[i] = (2*i+1)%23 ;
+
+  auto convert_to_vectorized =
+    [](const Vector<float> &in
+       , AlignedVector<VectorizedArray<float> > &out)
+  {
+    std::transform (in.begin(), in.end(), out.begin(),
+                    [](const auto &val)
+    {
+      return make_vectorized_array(val);
+    }) ;
+  } ;
+  AlignedVector<VectorizedArray<float> > v1(w1.size()), v2(w1.size()), v3(w1.size());
+  convert_to_vectorized (w1, v1) ;
+
+  constexpr unsigned int macro_size = VectorizedArray<float>::n_array_elements ;
+  Vector<float> vec_flat(v1.size()*macro_size) ;
+  std::array<unsigned int,macro_size> offsets ;
+  for (unsigned int i=0; i<macro_size; ++i)
+    offsets[i] = v1.size() * i ;
+  auto subtract_and_assign =
+    [](AlignedVector<VectorizedArray<float> > &lhs
+       , const AlignedVector<VectorizedArray<float> > &rhs)
+  {
+    std::transform (lhs.begin(), lhs.end(), rhs.begin(), lhs.begin(),
+                    [](const auto lval, const auto rval)
+    {
+      return lval - rval;
+    }) ;
+  } ;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  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();
+  deallog << "Verification of vmult and inverse: "
+          << (norm < 1e-3 ? 0. : norm)  << std::endl;
+
+  FullMatrix<float> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk)[0];
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0];
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii)[0] * mass[1](j,jj)[0]*mass[0](k,kk)[0]
+                      + mass[2](i,ii)[0] * (laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0]);
+  full.vmult(w2, w1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  norm = vec_flat.linfty_norm();
+  deallog << "Verification of vmult: " << (norm < 1e-4 ? 0. : norm) << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(w2, w1);
+  mat.apply_inverse(v2, v1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  norm = vec_flat.linfty_norm();
+  deallog << "Verification of inverse: " << (norm < 5e-3 ? 0. : norm) << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1,1>();
+  do_test<1,2>();
+  do_test<1,5>();
+  do_test<2,1>();
+  do_test<2,2>();
+  do_test<2,5>();
+  do_test<2,11>();
+  do_test<3,1>();
+  do_test<3,2>();
+  do_test<3,3>();
+  do_test<3,7>();
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_vectorized_03.with_lapack=true.output b/tests/lac/tensor_product_matrix_vectorized_03.with_lapack=true.output
new file mode 100644 (file)
index 0000000..d9e36d9
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0.00000
+DEAL::Verification of vmult: 0.00000
+DEAL::Verification of inverse: 0.00000
diff --git a/tests/lac/tensor_product_matrix_vectorized_04.cc b/tests/lac/tensor_product_matrix_vectorized_04.cc
new file mode 100644 (file)
index 0000000..38097a6
--- /dev/null
@@ -0,0 +1,145 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Same as previous tests except that we initialize the 'TensorProductMatrix'
+// with the same mass and derivative matrix in each tensor direction,
+// respectively.
+// Note, all lanes compute the same.
+
+#include "../tests.h"
+#include "../testmatrix.h"
+#include <deal.II/lac/tensor_product_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/base/vectorization.h>
+
+template <int dim>
+void do_test(const unsigned int size)
+{
+  deallog << "Testing dim=" << dim << ", degree=" << size << std::endl;
+  Table<2,VectorizedArray<double> > init_mass(size, size);
+  Table<2,VectorizedArray<double> > init_laplace(size, size);
+  for (unsigned int i=0; i<size; ++i)
+    {
+      init_mass(i,i) = 2./3.;
+      if (i > 0)
+        init_mass(i,i-1) = 1./6.;
+      if (i<size-1)
+        init_mass(i,i+1) = 1./6.;
+      init_laplace(i,i) = 2.;
+      if (i > 0)
+        init_laplace(i,i-1) = -1.;
+      if (i < size-1)
+        init_laplace(i,i+1) = -1.;
+    }
+  std::array<Table<2,VectorizedArray<double> >, dim> mass ;
+  std::array<Table<2,VectorizedArray<double> >, dim> laplace ;
+  std::fill (mass.begin(), mass.end(), init_mass) ;
+  std::fill (laplace.begin(), laplace.end(), init_laplace) ;
+
+  TensorProductMatrixSymmetricSum<dim,VectorizedArray<double> > mat;
+  mat.reinit(init_mass, init_laplace); // !!!
+
+  Vector<double> w1(mat.m()), w2(mat.m()) ;
+  for (unsigned int i=0; i<w1.size(); ++i)
+    w1[i] = (2*i+1)%23 ;
+
+  auto convert_to_vectorized =
+    [](const Vector<double> &in
+       , AlignedVector<VectorizedArray<double> > &out)
+  {
+    std::transform (in.begin(), in.end(), out.begin(),
+                    [](const auto &val)
+    {
+      return make_vectorized_array(val);
+    }) ;
+  } ;
+  AlignedVector<VectorizedArray<double> > v1(w1.size()), v2(w1.size()), v3(w1.size());
+  convert_to_vectorized (w1, v1) ;
+
+  constexpr unsigned int macro_size = VectorizedArray<double>::n_array_elements ;
+  Vector<double> vec_flat(v1.size()*macro_size) ;
+  std::array<unsigned int,macro_size> offsets ;
+  for (unsigned int i=0; i<macro_size; ++i)
+    offsets[i] = v1.size() * i ;
+  auto subtract_and_assign =
+    [](AlignedVector<VectorizedArray<double> > &lhs
+       , const AlignedVector<VectorizedArray<double> > &rhs)
+  {
+    std::transform (lhs.begin(), lhs.end(), rhs.begin(), lhs.begin(),
+                    [](const auto lval, const auto rval)
+    {
+      return lval - rval;
+    }) ;
+  } ;
+
+  mat.vmult(v2, v1);
+  mat.apply_inverse(v3, v2);
+  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;
+
+  FullMatrix<double> full(v1.size(), v1.size());
+  full = 0. ;
+  for (unsigned int dir = 0; dir<dim; ++dir)
+    for (unsigned int i=0, c=0; i<(dim>2?size:1); ++i)
+      for (unsigned int j=0; j<(dim>1?size:1); ++j)
+        for (unsigned int k=0; k<size; ++k, ++c)
+          for (unsigned int ii=0, cc=0; ii<(dim>2?size:1); ++ii)
+            for (unsigned int jj=0; jj<(dim>1?size:1); ++jj)
+              for (unsigned int kk=0; kk<size; ++kk, ++cc)
+                if (dim == 1)
+                  full(c,cc) = laplace[0](k,kk)[0];
+                else if (dim==2)
+                  full(c,cc) = laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0];
+                else if (dim==3)
+                  full(c,cc)
+                    = laplace[2](i,ii)[0] * mass[1](j,jj)[0]*mass[0](k,kk)[0]
+                      + mass[2](i,ii)[0] * (laplace[1](j,jj)[0]*mass[0](k,kk)[0] + mass[1](j,jj)[0]*laplace[0](k,kk)[0]);
+  full.vmult(w2, w1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of vmult: " << vec_flat.linfty_norm() << std::endl;
+
+  full.gauss_jordan();
+  full.vmult(w2, w1);
+  mat.apply_inverse(v2, v1);
+  convert_to_vectorized (w2, v3) ;
+  subtract_and_assign (v3, v2) ;
+  vectorized_transpose_and_store (false, mat.m(), v3.begin(), offsets.data(), vec_flat.begin()) ;
+  deallog << "Verification of inverse: " << vec_flat.linfty_norm() << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+
+  do_test<1>(1);
+  do_test<1>(2);
+  do_test<1>(5);
+  do_test<2>(1);
+  do_test<2>(2);
+  do_test<2>(5);
+  do_test<2>(11);
+  do_test<3>(1);
+  do_test<3>(2);
+  do_test<3>(3);
+  do_test<3>(7);
+
+  return 0;
+}
diff --git a/tests/lac/tensor_product_matrix_vectorized_04.with_lapack=true.output b/tests/lac/tensor_product_matrix_vectorized_04.with_lapack=true.output
new file mode 100644 (file)
index 0000000..4ef63ad
--- /dev/null
@@ -0,0 +1,45 @@
+
+DEAL::Testing dim=1, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=1, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=5
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=2, degree=11
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=1
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=2
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=3
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0
+DEAL::Testing dim=3, degree=7
+DEAL::Verification of vmult and inverse: 0
+DEAL::Verification of vmult: 0
+DEAL::Verification of inverse: 0

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.