*/
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
*/
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.
*/
*/
std::array<Table<2,Number>,dim> eigenvectors;
+private:
/**
- * An array for temporary data.
- */
+ * An array for temporary data.
+ */
mutable AlignedVector<Number> tmp_array;
/**
* 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) ;
*/
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.
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) ;
};
/**
* 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) ;
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) ;
};
* 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()) ;
-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
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) ;
}
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()),
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)) ;
}
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)) ;
}
{
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) ;
::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());
}
::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
::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) ;
}
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) ;
}
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)
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()) ;
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)
, 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)) ;
}
+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
::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());
}
::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());
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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