From: Julius Witte Date: Wed, 13 Sep 2017 14:28:39 +0000 (+0200) Subject: Added tests for all kinds or reinit() and vectorized template specialization of Tenso... X-Git-Tag: v9.0.0-rc1~1021^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6f057ed2ed9c5debcb8c0af877df4cefe8f844b8;p=dealii.git Added tests for all kinds or reinit() and vectorized template specialization of TensorProductMatrixSymmetricSum. Removed fill_data() from base class as we can access protected data members directly. --- diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 8f09f9dcb6..99af77b03e 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -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,dim> &mass_matrix, - const std::array,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 - 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,dim> eigenvectors; +private: /** - * An array for temporary data. - */ + * An array for temporary data. + */ mutable AlignedVector tmp_array; /** @@ -208,28 +189,28 @@ public: * immediately calling the corresponding reinit(). */ TensorProductMatrixSymmetricSum (const std::array, dim> &mass_matrix, - const std::array, dim> &derivative_matrix) ; - + const std::array, dim> &derivative_matrix) ; + /** * Constructor that is equivalent to the first constructor and * immediately calling the corresponding reinit(). */ TensorProductMatrixSymmetricSum (const std::array,dim> &mass_matrix, const std::array,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 &mass_matrix, - const FullMatrix &derivative_matrix) ; + const FullMatrix &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,dim> &mass_matrix, const std::array,dim> &derivative_matrix) ; @@ -240,7 +221,7 @@ public: */ void reinit (const std::array,dim> &mass_matrix, const std::array,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 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 > &mass_matrix, - const Table<2,VectorizedArray > &derivative_matrix) ; + const Table<2,VectorizedArray > &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 >,dim> &mass_matrix, const std::array >,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 void reinit_impl (MatrixArray &&mass_matrix, - MatrixArray &&derivative_matrix) ; + MatrixArray &&derivative_matrix) ; }; @@ -378,12 +359,13 @@ namespace * possible) */ template - 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 -template -inline -void -TensorProductMatrixSymmetricSumBase -::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(mass_matrices) ; - this->derivative_matrix = std::forward(derivative_matrices) ; - this->eigenvalues = std::forward(eigenvalues) ; - this->eigenvectors = std::forward(eigenvectors) ; -} - - - template inline unsigned int @@ -600,9 +560,9 @@ template inline TensorProductMatrixSymmetricSum ::TensorProductMatrixSymmetricSum (const std::array, dim> &mass_matrix, - const std::array, dim> &derivative_matrix) + const std::array, dim> &derivative_matrix) { - reinit_impl (mass_matrix, derivative_matrix) ; + reinit (mass_matrix, derivative_matrix) ; } @@ -635,13 +595,13 @@ inline void TensorProductMatrixSymmetricSum ::reinit_impl (MatrixArray &&mass_matrices_, - MatrixArray &&derivative_matrices_) + MatrixArray &&derivative_matrices_) { auto &&mass_matrices = std::forward(mass_matrices_) ; auto &&derivative_matrices = std::forward(derivative_matrices_) ; - - std::array,dim> eigenvectors ; - std::array, 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(size) == mass_matrices[dir].n_rows()), @@ -650,19 +610,15 @@ TensorProductMatrixSymmetricSum 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 (&(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 - ::fill_data (std::forward(mass_matrices), std::forward(derivative_matrices), - std::move(eigenvalues), std::move(eigenvectors)) ; } @@ -690,9 +646,9 @@ TensorProductMatrixSymmetricSum std::array,dim> deriv_copy ; std::transform (mass_matrix.cbegin(), mass_matrix.cend(), mass_copy.begin(), - [] (const FullMatrix& m) ->Table<2,Number> {return m;}) ; + [] (const FullMatrix &m) ->Table<2,Number> {return m;}) ; std::transform (derivative_matrix.cbegin(), derivative_matrix.cend(), deriv_copy.begin(), - [] (const FullMatrix& m) ->Table<2,Number> {return m;}) ; + [] (const FullMatrix &m) ->Table<2,Number> {return m;}) ; reinit_impl (std::move(mass_copy), std::move(deriv_copy)) ; } @@ -708,7 +664,7 @@ TensorProductMatrixSymmetricSum { std::array,dim> mass_matrices ; std::array,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 ::vmult (Vector &dst, const Vector &src) const { - AssertDimension(dst.size(), Utilities::fixed_power(this->eigenvalues[0].size())); - AssertDimension(src.size(), Utilities::fixed_power(this->eigenvalues[0].size())); + AssertDimension(dst.size(), this->m()) ; + AssertDimension(src.size(), this->n()) ; TensorProductMatrixSymmetricSumBase::vmult (dst.begin(), src.begin()); } @@ -738,63 +694,14 @@ TensorProductMatrixSymmetricSum ::apply_inverse (Vector &dst, const Vector &src) const { - AssertDimension (dst.size(), Utilities::fixed_power(this->eigenvalues[0].size())); - AssertDimension (src.size(), Utilities::fixed_power(this->eigenvalues[0].size())); + AssertDimension (dst.size(), this->n()) ; + AssertDimension (src.size(), this->m()) ; TensorProductMatrixSymmetricSumBase::apply_inverse (dst.begin(), src.begin()); } -// template -// inline -// FullMatrix -// TensorProductMatrixSymmetricSum -// ::get_full_matrix () const -// { -// const auto& mass_matrix = TensorProductMatrixSymmetricSumBase::mass_matrix ; -// const auto& derivative_matrix = this->derivative_matrix ; -// const auto& eigenvalues = this->eigenvalues ; - -// FullMatrix matrix {Utilities::fixed_power(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 inline @@ -811,7 +718,7 @@ TensorProductMatrixSymmetricSum,size> ::TensorProductMatrixSymmetricSum (const std::array >,dim> &mass_matrix, const std::array >,dim> &derivative_matrix) { - reinit_impl (mass_matrix, derivative_matrix) ; + reinit (mass_matrix, derivative_matrix) ; } @@ -820,15 +727,9 @@ template inline TensorProductMatrixSymmetricSum,size> ::TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray > &mass_matrix, - const Table<2,VectorizedArray > &derivative_matrix) + const Table<2,VectorizedArray > &derivative_matrix) { - std::array >,dim> mass_matrices ; - std::array >,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,size> ::reinit_impl (MatrixArray &&mass_matrices_, - MatrixArray &&derivative_matrices_) + MatrixArray &&derivative_matrices_) { auto &&mass_matrix = std::forward(mass_matrices_) ; auto &&derivative_matrix = std::forward(derivative_matrices_) ; - std::array >,dim> eigenvectors ; - std::array >, dim> eigenvalues ; - + this->mass_matrix = mass_matrix ; + this->derivative_matrix = derivative_matrix ; + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + const unsigned int nm_flat_size + = (size > 0) + ? (Utilities::fixed_int_power::value + * Utilities::fixed_int_power::value * macro_size) + : (Utilities::fixed_power(mass_matrix[0].n_rows()) + * Utilities::fixed_power(mass_matrix[0].n_rows()) * macro_size) ; + const unsigned int n_flat_size + = (size > 0) + ? Utilities::fixed_int_power::value * macro_size + : Utilities::fixed_power(mass_matrix[0].n_rows()) * macro_size ; std::vector mass_matrix_flat ; std::vector deriv_matrix_flat ; std::vector eigenvalues_flat ; std::vector 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 offsets_nm ; std::array offsets_n ; for (int dir = 0; dir < dim; ++dir) @@ -866,14 +781,15 @@ TensorProductMatrixSymmetricSum,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,size> const Number *deriv_cbegin = deriv_matrix_flat.data() ; Number *eigenvec_begin = eigenvectors_flat.data() ; Number *eigenval_begin = eigenvalues_flat.data() ; + spectral_assembly (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,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,size> - ::fill_data (std::forward(mass_matrix), std::forward(derivative_matrix), - std::move(eigenvalues), std::move(eigenvectors)) ; } @@ -923,6 +836,24 @@ TensorProductMatrixSymmetricSum,size> +template +inline +void +TensorProductMatrixSymmetricSum,size> +::reinit (const Table<2,VectorizedArray > &mass_matrix, + const Table<2,VectorizedArray > &derivative_matrix) +{ + std::array >,dim> mass_matrices ; + std::array >,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 inline void @@ -930,8 +861,8 @@ TensorProductMatrixSymmetricSum,size> ::vmult (AlignedVector > &dst, const AlignedVector > &src) const { - AssertDimension(dst.size(), Utilities::fixed_power(this->eigenvalues[0].size())); - AssertDimension(src.size(), Utilities::fixed_power(this->eigenvalues[0].size())); + AssertDimension(dst.size(), this->m()) ; + AssertDimension(src.size(), this->n()) ; TensorProductMatrixSymmetricSumBase,size>::vmult (dst.begin(), src.begin()); } @@ -944,8 +875,8 @@ TensorProductMatrixSymmetricSum,size> ::apply_inverse (AlignedVector > &dst, const AlignedVector > &src) const { - AssertDimension (dst.size(), Utilities::fixed_power(this->eigenvalues[0].size())); - AssertDimension (src.size(), Utilities::fixed_power(this->eigenvalues[0].size())); + AssertDimension (dst.size(), this->n()) ; + AssertDimension (src.size(), this->m()) ; TensorProductMatrixSymmetricSumBase,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 index 0000000000..4fc47380b7 --- /dev/null +++ b/tests/lac/tensor_product_matrix_04.cc @@ -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 +#include +#include + + +template +void do_test(const unsigned int size) +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + FullMatrix init_mass(size, size); + FullMatrix init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6. ; + if (i 0) + init_laplace(i,i-1) = -1. ; + if (i < size-1) + init_laplace(i,i+1) = -1. ; + } + + std::array, dim> mass ; + std::array, dim> laplace ; + for (unsigned int dir = 0; dir mat; + mat.reinit(mass, laplace); + Vector v1(mat.m()), v2(mat.m()), v3(mat.m()); + for (unsigned int i=0; i full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(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 index 0000000000..bde376ccf4 --- /dev/null +++ b/tests/lac/tensor_product_matrix_04.with_lapack=true.output @@ -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 index 0000000000..d0f84f69cd --- /dev/null +++ b/tests/lac/tensor_product_matrix_05.cc @@ -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 +#include +#include + + +template +void do_test() +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + FullMatrix init_mass(size, size); + FullMatrix init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6. ; + if (i 0) + init_laplace(i,i-1) = -1. ; + if (i < size-1) + init_laplace(i,i+1) = -1. ; + } + + std::array, dim> mass ; + std::array, dim> laplace ; + for (unsigned int dir = 0; dir mat; + mat.reinit(mass, laplace); + Vector v1(mat.m()), v2(mat.m()), v3(mat.m()); + for (unsigned int i=0; i full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(); + 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 index 0000000000..bde376ccf4 --- /dev/null +++ b/tests/lac/tensor_product_matrix_05.with_lapack=true.output @@ -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 index 0000000000..2ec0b53419 --- /dev/null +++ b/tests/lac/tensor_product_matrix_06.cc @@ -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 +#include +#include + + +template +void do_test() +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + FullMatrix init_mass(size, size); + FullMatrix init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6. ; + if (i 0) + init_laplace(i,i-1) = -1. ; + if (i < size-1) + init_laplace(i,i+1) = -1. ; + } + + std::array, dim> mass ; + std::array, dim> laplace ; + for (unsigned int dir = 0; dir mat; + mat.reinit(mass, laplace); + Vector v1(mat.m()), v2(mat.m()), v3(mat.m()); + for (unsigned int i=0; i full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(); + 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 index 0000000000..2fc78058b3 --- /dev/null +++ b/tests/lac/tensor_product_matrix_06.with_lapack=true.output @@ -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 index 0000000000..ec37cb87c2 --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_01.cc @@ -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'. +// Note, all lanes compute the same. + +#include "../tests.h" +#include "../testmatrix.h" +#include +#include +#include +#include + +template +void do_test(const unsigned int size) +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + Table<2,VectorizedArray > init_mass(size, size); + Table<2,VectorizedArray > init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6.; + if (i 0) + init_laplace(i,i-1) = -1.; + if (i < size-1) + init_laplace(i,i+1) = -1.; + } + std::array >, dim> mass ; + std::array >, dim> laplace ; + for (unsigned int dir = 0; dir(4./3.) ; + init_laplace(i,i) *= make_vectorized_array(5./4.) ; + } + mass[dir] = init_mass ; + laplace[dir] = init_laplace ; + } + TensorProductMatrixSymmetricSum > mat; + mat.reinit(mass, laplace); + + Vector w1(mat.m()), w2(mat.m()) ; + for (unsigned int i=0; i &in + , AlignedVector > &out) + { + std::transform (in.begin(), in.end(), out.begin(), + [](const auto &val) + { + return make_vectorized_array(val); + }) ; + } ; + AlignedVector > v1(w1.size()), v2(w1.size()), v3(w1.size()); + convert_to_vectorized (w1, v1) ; + + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + Vector vec_flat(v1.size()*macro_size) ; + std::array offsets ; + for (unsigned int i=0; i > &lhs + , const AlignedVector > &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 full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(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 index 0000000000..4ef63ad3dc --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_01.with_lapack=true.output @@ -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 index 0000000000..db5976138f --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_02.cc @@ -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'. +// Note, all lanes compute the same. + +#include "../tests.h" +#include "../testmatrix.h" +#include +#include +#include +#include + +template +void do_test() +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + Table<2,VectorizedArray > init_mass(size, size); + Table<2,VectorizedArray > init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6.; + if (i 0) + init_laplace(i,i-1) = -1.; + if (i < size-1) + init_laplace(i,i+1) = -1.; + } + std::array >, dim> mass ; + std::array >, dim> laplace ; + for (unsigned int dir = 0; dir(4./3.) ; + init_laplace(i,i) *= make_vectorized_array(5./4.) ; + } + mass[dir] = init_mass ; + laplace[dir] = init_laplace ; + } + TensorProductMatrixSymmetricSum > mat; + mat.reinit(mass, laplace); + + Vector w1(mat.m()), w2(mat.m()) ; + for (unsigned int i=0; i &in + , AlignedVector > &out) + { + std::transform (in.begin(), in.end(), out.begin(), + [](const auto &val) + { + return make_vectorized_array(val); + }) ; + } ; + AlignedVector > v1(w1.size()), v2(w1.size()), v3(w1.size()); + convert_to_vectorized (w1, v1) ; + + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + Vector vec_flat(v1.size()*macro_size) ; + std::array offsets ; + for (unsigned int i=0; i > &lhs + , const AlignedVector > &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 full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(); + 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 index 0000000000..4ef63ad3dc --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_02.with_lapack=true.output @@ -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 index 0000000000..0e73c3aa32 --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_03.cc @@ -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'. +// Note, all lanes compute the same. + +#include "../tests.h" +#include "../testmatrix.h" +#include +#include +#include +#include + +template +void do_test() +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + Table<2,VectorizedArray > init_mass(size, size); + Table<2,VectorizedArray > init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6.; + if (i 0) + init_laplace(i,i-1) = -1.; + if (i < size-1) + init_laplace(i,i+1) = -1.; + } + std::array >, dim> mass ; + std::array >, dim> laplace ; + for (unsigned int dir = 0; dir(4./3.) ; + init_laplace(i,i) *= make_vectorized_array(5./4.) ; + } + mass[dir] = init_mass ; + laplace[dir] = init_laplace ; + } + TensorProductMatrixSymmetricSum > mat; + mat.reinit(mass, laplace); + + Vector w1(mat.m()), w2(mat.m()) ; + for (unsigned int i=0; i &in + , AlignedVector > &out) + { + std::transform (in.begin(), in.end(), out.begin(), + [](const auto &val) + { + return make_vectorized_array(val); + }) ; + } ; + AlignedVector > v1(w1.size()), v2(w1.size()), v3(w1.size()); + convert_to_vectorized (w1, v1) ; + + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + Vector vec_flat(v1.size()*macro_size) ; + std::array offsets ; + for (unsigned int i=0; i > &lhs + , const AlignedVector > &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 full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(); + 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 index 0000000000..d9e36d9e7a --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_03.with_lapack=true.output @@ -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 index 0000000000..38097a6c65 --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_04.cc @@ -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 +#include +#include +#include + +template +void do_test(const unsigned int size) +{ + deallog << "Testing dim=" << dim << ", degree=" << size << std::endl; + Table<2,VectorizedArray > init_mass(size, size); + Table<2,VectorizedArray > init_laplace(size, size); + for (unsigned int i=0; i 0) + init_mass(i,i-1) = 1./6.; + if (i 0) + init_laplace(i,i-1) = -1.; + if (i < size-1) + init_laplace(i,i+1) = -1.; + } + std::array >, dim> mass ; + std::array >, dim> laplace ; + std::fill (mass.begin(), mass.end(), init_mass) ; + std::fill (laplace.begin(), laplace.end(), init_laplace) ; + + TensorProductMatrixSymmetricSum > mat; + mat.reinit(init_mass, init_laplace); // !!! + + Vector w1(mat.m()), w2(mat.m()) ; + for (unsigned int i=0; i &in + , AlignedVector > &out) + { + std::transform (in.begin(), in.end(), out.begin(), + [](const auto &val) + { + return make_vectorized_array(val); + }) ; + } ; + AlignedVector > v1(w1.size()), v2(w1.size()), v3(w1.size()); + convert_to_vectorized (w1, v1) ; + + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + Vector vec_flat(v1.size()*macro_size) ; + std::array offsets ; + for (unsigned int i=0; i > &lhs + , const AlignedVector > &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 full(v1.size(), v1.size()); + full = 0. ; + for (unsigned int dir = 0; dir2?size:1); ++i) + for (unsigned int j=0; j<(dim>1?size:1); ++j) + for (unsigned int k=0; k2?size:1); ++ii) + for (unsigned int jj=0; jj<(dim>1?size:1); ++jj) + for (unsigned int kk=0; kk(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 index 0000000000..4ef63ad3dc --- /dev/null +++ b/tests/lac/tensor_product_matrix_vectorized_04.with_lapack=true.output @@ -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