From: Julius Witte Date: Tue, 5 Sep 2017 15:18:28 +0000 (+0200) Subject: Split up of TensorProductMatrixSymmetricSum into base and derived class. Added templa... X-Git-Tag: v9.0.0-rc1~1021^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=172e35f989f3bc1aa8f9d7e5c8ac38dba3111425;p=dealii.git Split up of TensorProductMatrixSymmetricSum into base and derived class. Added template specialization for arithmetic template Number being VectorizedArray. --- diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index b4804e1054..8f09f9dcb6 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -97,20 +97,51 @@ template class FullMatrix; * @author Martin Kronbichler, 2017 */ template -class TensorProductMatrixSymmetricSum +class TensorProductMatrixSymmetricSumBase { public: + /** + * Returns the number of rows of this matrix, given by the dim-th power of + * the size of the 1D matrices passed to the constructor. + */ + unsigned int m () const; + + /** + * Returns the number of columns of this matrix, given by the dim-th power + * of the size of the 1D matrices passed to the constructor. + */ + unsigned int n () const; + + /** + * Implements a matrix-vector product with the underlying matrix as + * described in the main documentation of this class. Same as the other + * vmult() function, but operating on plain pointers rather than a vector + * (no check of array bounds possible). + */ + 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 + * apply_inverse() function, but operating on plain pointers rather than a + * vector (no check of array bounds possible). + */ + void apply_inverse (Number *dst, + const Number *src) const; + +protected: /** * Constructor. */ - TensorProductMatrixSymmetricSum() = default; + TensorProductMatrixSymmetricSumBase () = default ; /** * Constructor that is equivalent to the previous constructor and * immediately calling reinit(). */ - TensorProductMatrixSymmetricSum(const FullMatrix &mass_matrix, - const FullMatrix &derivative_matrix); + 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 @@ -118,36 +149,114 @@ public: * and positive definite and $A$ to be symmetric and invertible but not * necessarily positive defininte. */ - void reinit (const FullMatrix &mass_matrix, - const FullMatrix &derivative_matrix); + template + void + fill_data (MatrixArray&& mass_matrices, + MatrixArray&& derivative_matrices, + EigenvalueType&& eigenvalues, + EigenvectorType&& eigenvectors) ; /** - * Returns the number of rows of this matrix, given by the dim-th power of - * the size of the 1D matrices passed to the constructor. + * A copy of the @p mass_matrix object passed to the reinit() method. */ - unsigned int m() const; + std::array,dim> mass_matrix; /** - * Returns the number of columns of this matrix, given by the dim-th power - * of the size of the 1D matrices passed to the constructor. + * A copy of the @p derivative_matrix object passed to the reinit() method. */ - unsigned int n() const; + std::array,dim> derivative_matrix; /** - * Implements a matrix-vector product with the underlying matrix as - * described in the main documentation of this class. + * A vector containing the generalized eigenvalues of A s = lambda B s. */ - void vmult (Vector &dst, - const Vector &src) const; + std::array,dim> eigenvalues; + + /** + * The matrix containing the generalized eigenvectors. + */ + std::array,dim> eigenvectors; + + /** + * An array for temporary data. + */ + mutable AlignedVector tmp_array; + + /** + * A mutex that guards access to the array @p tmp_array. + */ + mutable Threads::Mutex mutex; +}; + + + +/** + * ... new TensorProductMatrixSymmetricSum using the base class as tensor product + * container and interface to arithmetic operations for a generic Number type ... + */ +template +class TensorProductMatrixSymmetricSum + : public TensorProductMatrixSymmetricSumBase +{ +public: + /** + * Constructor. + */ + TensorProductMatrixSymmetricSum () ; + + /** + * Constructor that is equivalent to the previous 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(). + */ + 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(). + */ + TensorProductMatrixSymmetricSum (const FullMatrix &mass_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. + */ + void reinit (const std::array,dim> &mass_matrix, + const std::array,dim> &derivative_matrix) ; + + /** + * Equivalent to the previous reinit() unless that the mass and derivative + * matrices are passed by Table instead of FullMatrix. + */ + 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. + * 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. + */ + void reinit (const FullMatrix &mass_matrix, + const FullMatrix &derivative_matrix) ; /** * Implements a matrix-vector product with the underlying matrix as - * described in the main documentation of this class. Same as the other - * vmult() function, but operating on plain pointers rather than a vector - * (no check of array bounds possible). + * described in the main documentation of this class. */ - void vmult (Number *dst, - const Number *src) const; + void vmult (Vector &dst, + const Vector &src) const; /** * Implements a matrix-vector product with the underlying matrix as @@ -157,44 +266,100 @@ public: const Vector &src) const; /** - * Implements a matrix-vector product with the underlying matrix as - * described in the main documentation of this class. Same as the other - * apply_inverse() function, but operating on plain pointers rather than a - * vector (no check of array bounds possible). + * ... for compability to MappingQGeneric */ - void apply_inverse (Number *dst, - const Number *src) const; + using TensorProductMatrixSymmetricSumBase::vmult ; + + /** + * ... for compability to MappingQGeneric + */ + using TensorProductMatrixSymmetricSumBase::apply_inverse ; private: /** - * A copy of the @p mass_matrix object passed to the reinit() method. + * 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. */ - FullMatrix mass_matrix; + template + void reinit_impl (MatrixArray &&mass_matrix, + MatrixArray &&derivative_matrix) ; +}; + +/** + * ... same as previous class but based on a vectorized value type, namely + * VectorizedArray ... + */ +template +class TensorProductMatrixSymmetricSum,size> + : public TensorProductMatrixSymmetricSumBase,size> +{ +public: /** - * A copy of the @p derivative_matrix object passed to the reinit() method. + * Constructor. */ - FullMatrix derivative_matrix; + TensorProductMatrixSymmetricSum () ; /** - * A vector containing the generalized eigenvalues of A s = lambda B s. + * Constructor that is equivalent to the previous constructor and + * immediately calling reinit(). */ - AlignedVector eigenvalues; + TensorProductMatrixSymmetricSum (const std::array >,dim> &mass_matrix, + const std::array >,dim> &derivative_matrix) ; /** - * The matrix containing the generalized eigenvectors. + * Constructor that is equivalent to the first constructor and + * immediately calling the corresponding reinit(). */ - Table<2,Number> eigenvectors; + TensorProductMatrixSymmetricSum (const Table<2,VectorizedArray > &mass_matrix, + const Table<2,VectorizedArray > &derivative_matrix) ; /** - * An array for temporary data. + * 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. */ - mutable AlignedVector tmp_array; + void reinit (const std::array >,dim> &mass_matrix, + const std::array >,dim> &derivative_matrix) ; /** - * A mutex that guards access to the array @p tmp_array. + * Initializes the same mass matrix $M$ and derivative matrix $A$ to the given array + * of mass matrices and array of derivative matrices, respectively. + * 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. */ - mutable Threads::Mutex mutex; + void reinit (const Table<2,VectorizedArray > &mass_matrix, + const Table<2,VectorizedArray > &derivative_matrix) ; + + /** + * Implements a matrix-vector product with the underlying matrix as + * described in the main documentation of this class. + */ + void vmult (AlignedVector > &dst, + const AlignedVector > &src) const ; + + /** + * Implements a matrix-vector product with the underlying matrix as + * described in the main documentation of this class. + */ + void apply_inverse (AlignedVector > &dst, + const AlignedVector > &src) const ; + +private: + /** + * 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) ; }; @@ -202,51 +367,72 @@ private: #ifndef DOXYGEN - -template -inline -TensorProductMatrixSymmetricSum -::TensorProductMatrixSymmetricSum(const FullMatrix &mass_matrix, - const FullMatrix &derivative_matrix) +namespace { - reinit(mass_matrix, derivative_matrix); + /** + * Compute generalized eigenvalues and eigenvectors of the real + * generalized symmetric eigenproblem $M v = \lambda A v$. Since we are + * operating on plain pointers we require the size of the matrices beforehand. + * Note that the data arrays for the eigenvalues and eigenvectors + * have to be initialized to a proper size, too. (no check of array bounds + * possible) + */ + template + 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()) ; + + auto &&transpose_fill_nm + = [](Number *out, const Number *in, const unsigned int n, const unsigned int m) + { + for (unsigned int mm = 0; mm < m; ++mm) + for (unsigned int nn = 0; nn < n; ++nn) + out[mm+nn*m] = *(in++) ; + }; + + std::vector > eigenvecs(n_rows) ; + LAPACKFullMatrix mass_copy(n_rows, n_cols) ; + LAPACKFullMatrix deriv_copy(n_rows, n_cols) ; + + transpose_fill_nm (&(mass_copy(0,0)), mass_matrix, n_rows, n_cols) ; + transpose_fill_nm (&(deriv_copy(0,0)), derivative_matrix, n_rows, n_cols) ; + + deriv_copy.compute_generalized_eigenvalues_symmetric (mass_copy, eigenvecs); + AssertDimension (eigenvecs.size(), n_rows) ; + for (unsigned int i=0; i +template inline void -TensorProductMatrixSymmetricSum -::reinit(const FullMatrix &mass_matrix, - const FullMatrix &derivative_matrix) +TensorProductMatrixSymmetricSumBase +::fill_data (MatrixArray&& mass_matrices, + MatrixArray&& derivative_matrices, + EigenvalueType&& eigenvalues, + EigenvectorType&& eigenvectors) { - Assert(size == -1 || - (size > 0 && static_cast(size) == mass_matrix.m()), - ExcDimensionMismatch(size, mass_matrix.m())); - AssertDimension(mass_matrix.m(), mass_matrix.n()); - AssertDimension(mass_matrix.m(), derivative_matrix.m()); - AssertDimension(mass_matrix.m(), derivative_matrix.n()); - - this->mass_matrix = mass_matrix; - this->derivative_matrix = derivative_matrix; - - std::vector > eigenvecs(mass_matrix.m()); - LAPACKFullMatrix mass_copy(mass_matrix.m(), mass_matrix.n()); - LAPACKFullMatrix deriv_copy(derivative_matrix.m(), derivative_matrix.n()); - mass_copy = mass_matrix; - deriv_copy = derivative_matrix; - - deriv_copy.compute_generalized_eigenvalues_symmetric(mass_copy, eigenvecs); - AssertDimension(eigenvecs.size(), mass_matrix.m()); - eigenvectors.reinit(mass_matrix.m(), mass_matrix.m()); - for (unsigned int i=0; imass_matrix = std::forward(mass_matrices) ; + this->derivative_matrix = std::forward(derivative_matrices) ; + this->eigenvalues = std::forward(eigenvalues) ; + this->eigenvectors = std::forward(eigenvectors) ; } @@ -254,9 +440,12 @@ TensorProductMatrixSymmetricSum template inline unsigned int -TensorProductMatrixSymmetricSum::m() const +TensorProductMatrixSymmetricSumBase::m() const { - return Utilities::fixed_power(mass_matrix.m()); + unsigned int m = mass_matrix[0].n_rows() ; + for (unsigned int d = 1; d < dim; ++d) + m *= mass_matrix[d].n_rows() ; + return m ; } @@ -264,9 +453,12 @@ TensorProductMatrixSymmetricSum::m() const template inline unsigned int -TensorProductMatrixSymmetricSum::n() const +TensorProductMatrixSymmetricSumBase::n() const { - return Utilities::fixed_power(mass_matrix.n()); + unsigned int n = mass_matrix[0].n_cols() ; + for (unsigned int d = 1; d < dim; ++d) + n *= mass_matrix[d].n_cols() ; + return n ; } @@ -274,67 +466,54 @@ TensorProductMatrixSymmetricSum::n() const template inline void -TensorProductMatrixSymmetricSum -::vmult(Vector &dst, - const Vector &src) const -{ - AssertDimension(dst.size(), Utilities::fixed_power(eigenvalues.size())); - AssertDimension(src.size(), Utilities::fixed_power(eigenvalues.size())); - vmult(dst.begin(), src.begin()); -} - - - -template -inline -void -TensorProductMatrixSymmetricSum -::apply_inverse(Vector &dst, - const Vector &src) const -{ - AssertDimension(dst.size(), Utilities::fixed_power(eigenvalues.size())); - AssertDimension(src.size(), Utilities::fixed_power(eigenvalues.size())); - apply_inverse(dst.begin(), src.begin()); -} - - - -template -inline -void -TensorProductMatrixSymmetricSum +TensorProductMatrixSymmetricSumBase ::vmult(Number *dst, const Number *src) const { Threads::Mutex::ScopedLock lock(this->mutex); - const unsigned int n = Utilities::fixed_power(size > 0 ? size : eigenvalues.size()); + const unsigned int n = Utilities::fixed_power(size > 0 ? size : eigenvalues[0].size()); tmp_array.resize_fast(n*2); - const int kernel_size = size > 0 ? size-1 : -1; + constexpr int kernel_size = size > 0 ? size-1 : -1; internal::EvaluatorTensorProduct - eval(AlignedVector(), AlignedVector(), - AlignedVector(), mass_matrix.m()-1, mass_matrix.m()); - const Number *A = &derivative_matrix(0,0); - const Number *M = &mass_matrix(0,0); + eval(AlignedVector {}, AlignedVector {}, + AlignedVector {}, mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows()); Number *t = tmp_array.begin(); + if (dim == 1) - eval.template apply<0, true, false>(A, src, dst); + { + const Number *A = &derivative_matrix[0](0,0); + eval.template apply<0, false, false> (A, src, dst); + } + else if (dim == 2) { - eval.template apply<0, true, false>(M, src, t); - eval.template apply<1, true, false>(A, t, dst); - eval.template apply<0, true, false>(A, src, t); - eval.template apply<1, true, true> (M, t, dst); + const Number *A0 = &derivative_matrix[0](0,0); + const Number *M0 = &mass_matrix[0](0,0); + const Number *A1 = &derivative_matrix[1](0,0); + const Number *M1 = &mass_matrix[1](0,0); + eval.template apply<0, false, false> (M0, src, t); + eval.template apply<1, false, false> (A1, t, dst); + eval.template apply<0, false, false> (A0, src, t); + eval.template apply<1, false, true> (M1, t, dst); } + else if (dim == 3) { - eval.template apply<0, true, false>(M, src, t+n); - eval.template apply<1, true, false>(M, t+n, t); - eval.template apply<2, true, false>(A, t, dst); - eval.template apply<1, true, false>(A, t+n, t); - eval.template apply<0, true, false>(A, src, t+n); - eval.template apply<1, true, true> (M, t+n, t); - eval.template apply<2, true, true> (M, t, dst); + const Number *A0 = &derivative_matrix[0](0,0); + const Number *M0 = &mass_matrix[0](0,0); + const Number *A1 = &derivative_matrix[1](0,0); + const Number *M1 = &mass_matrix[1](0,0); + const Number *A2 = &derivative_matrix[2](0,0); + const Number *M2 = &mass_matrix[2](0,0); + eval.template apply<0, false, false> (M0, src, t+n); + eval.template apply<1, false, false> (M1, t+n, t); + eval.template apply<2, false, false> (A2, t, dst); + eval.template apply<1, false, false> (A1, t+n, t); + eval.template apply<0, false, false> (A0, src, t+n); + eval.template apply<1, false, true> (M1, t+n, t); + eval.template apply<2, false, true> (M2, t, dst); } + else AssertThrow(false, ExcNotImplemented()); } @@ -344,55 +523,430 @@ TensorProductMatrixSymmetricSum template inline void -TensorProductMatrixSymmetricSum +TensorProductMatrixSymmetricSumBase ::apply_inverse(Number *dst, const Number *src) const { Threads::Mutex::ScopedLock lock(this->mutex); - const unsigned int n = size > 0 ? size : eigenvalues.size(); - tmp_array.resize_fast(Utilities::fixed_power(n)); - const int kernel_size = size > 0 ? size-1 : -1; + const unsigned int n = size > 0 ? size : eigenvalues[0].size(); + tmp_array.resize_fast (Utilities::fixed_power(n)); + constexpr int kernel_size = size > 0 ? size-1 : -1; internal::EvaluatorTensorProduct eval(AlignedVector(), AlignedVector(), - AlignedVector(), mass_matrix.m()-1, mass_matrix.m()); - const Number *S = &eigenvectors(0,0); + AlignedVector(), mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows()); Number *t = tmp_array.begin(); - switch (dim) + // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index' + // --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src, + // --> apply<.,false,.> (S,src,dst) calculates dst = S * src, + // while the eigenvectors are stored column-wise in S, i.e. + // rows correspond to dofs whereas columns to eigenvalue indices! + if (dim == 1) { - case 1: + const Number *S = &eigenvectors[0](0,0); eval.template apply<0, true, false> (S, src, t); for (unsigned int i=0; i (S, t, dst); - break; + } + + else if (dim == 2) + { + const Number *S0 = &(eigenvectors[0](0,0)); + const Number *S1 = &(eigenvectors[1](0,0)); + eval.template apply<0, true, false> (S0, src, t); + eval.template apply<1, true, false> (S1, t, dst); + for (unsigned int i1=0, c=0; i1 (S0, dst, t); + eval.template apply<1, false, false> (S1, t, dst); + } + + else if (dim == 3) + { + const Number *S0 = &eigenvectors[0](0,0); + const Number *S1 = &eigenvectors[1](0,0); + const Number *S2 = &eigenvectors[2](0,0); + eval.template apply<0, true, false> (S0, src, t); + eval.template apply<1, true, false> (S1, t, dst); + eval.template apply<2, true, false> (S2, dst, t); + for (unsigned int i2=0, c=0; i2 (S0, t, dst); + eval.template apply<1, false, false> (S1, dst, t); + eval.template apply<2, false, false> (S2, t, dst); + } + + else + Assert(false, ExcNotImplemented()); +} + + +// ------------------------------ TensorProductMatrixSymmetricSum ------------------------------ + +template +inline +TensorProductMatrixSymmetricSum +::TensorProductMatrixSymmetricSum () + : TensorProductMatrixSymmetricSumBase() +{} + + + +template +inline +TensorProductMatrixSymmetricSum +::TensorProductMatrixSymmetricSum (const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) +{ + reinit_impl (mass_matrix, derivative_matrix) ; +} + + + +template +inline +TensorProductMatrixSymmetricSum +::TensorProductMatrixSymmetricSum(const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) +{ + reinit (mass_matrix, derivative_matrix) ; +} + + + +template +inline +TensorProductMatrixSymmetricSum +::TensorProductMatrixSymmetricSum (const FullMatrix &mass_matrix, + const FullMatrix &derivative_matrix) +{ + reinit (mass_matrix, derivative_matrix) ; +} + + + +template +template +inline +void +TensorProductMatrixSymmetricSum +::reinit_impl (MatrixArray &&mass_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 ; + for (int dir = 0; dir < dim; ++dir) + { + Assert (size == -1 || (size > 0 && static_cast(size) == mass_matrices[dir].n_rows()), + ExcDimensionMismatch(size, mass_matrices[dir].n_rows())); + AssertDimension (mass_matrices[dir].n_rows(), mass_matrices[dir].n_cols()); + 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()) ; + 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))) ; + } + + TensorProductMatrixSymmetricSumBase + ::fill_data (std::forward(mass_matrices), std::forward(derivative_matrices), + std::move(eigenvalues), std::move(eigenvectors)) ; +} + + + +template +inline +void +TensorProductMatrixSymmetricSum +::reinit (const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) +{ + reinit_impl (mass_matrix, derivative_matrix) ; +} + + + +template +inline +void +TensorProductMatrixSymmetricSum +::reinit (const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) +{ + std::array,dim> mass_copy ; + std::array,dim> deriv_copy ; + + std::transform (mass_matrix.cbegin(), mass_matrix.cend(), mass_copy.begin(), + [] (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;}) ; + + reinit_impl (std::move(mass_copy), std::move(deriv_copy)) ; +} + + + +template +inline +void +TensorProductMatrixSymmetricSum +::reinit (const FullMatrix &mass_matrix, + const FullMatrix &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 +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())); + TensorProductMatrixSymmetricSumBase::vmult (dst.begin(), src.begin()); +} + + + +template +inline +void +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())); + 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 ------------------------------ + +template +inline +TensorProductMatrixSymmetricSum,size> +::TensorProductMatrixSymmetricSum () + : TensorProductMatrixSymmetricSumBase,size>() +{} + + + +template +inline +TensorProductMatrixSymmetricSum,size> +::TensorProductMatrixSymmetricSum (const std::array >,dim> &mass_matrix, + const std::array >,dim> &derivative_matrix) +{ + reinit_impl (mass_matrix, derivative_matrix) ; +} + + + +template +inline +TensorProductMatrixSymmetricSum,size> +::TensorProductMatrixSymmetricSum (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 (mass_matrices, derivative_matrices) ; +} - case 2: - eval.template apply<0, true, false> (S, src, t); - eval.template apply<1, true, false> (S, t, dst); - for (unsigned int i=0, c=0; i (S, dst, t); - eval.template apply<0, false, false> (S, t, dst); - break; - case 3: - eval.template apply<0, true, false> (S, src, t); - eval.template apply<1, true, false> (S, t, dst); - eval.template apply<2, true, false> (S, dst, t); - for (unsigned int i=0, c=0; i (S, t, dst); - eval.template apply<1, false, false> (S, dst, t); - eval.template apply<0, false, false> (S, t, dst); - break; - default: - Assert(false, ExcNotImplemented()); +template +template +inline +void +TensorProductMatrixSymmetricSum,size> +::reinit_impl (MatrixArray &&mass_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 ; + + constexpr unsigned int macro_size = VectorizedArray::n_array_elements ; + + std::vector mass_matrix_flat ; + std::vector deriv_matrix_flat ; + std::vector eigenvalues_flat ; + std::vector eigenvectors_flat ; + std::array offsets_nm ; + std::array offsets_n ; + for (int dir = 0; dir < dim; ++dir) + { + Assert (size == -1 || + (size > 0 && static_cast(size) == mass_matrix[dir].n_rows()), + ExcDimensionMismatch(size, mass_matrix[dir].n_rows())); + AssertDimension (mass_matrix[dir].n_rows(), mass_matrix[dir].n_cols()); + AssertDimension (mass_matrix[dir].n_rows(), derivative_matrix[dir].n_rows()); + AssertDimension (mass_matrix[dir].n_rows(), derivative_matrix[dir].n_cols()); + + 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++);}) ; + + vectorized_transpose_and_store (false, nm, &(mass_matrix[dir](0,0)) + , offsets_nm.cbegin(), mass_matrix_flat.data()) ; + vectorized_transpose_and_store (false, nm, &(derivative_matrix[dir](0,0)) + , offsets_nm.cbegin(), deriv_matrix_flat.data()) ; + + const Number *mass_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 (mass_cbegin, deriv_cbegin, n_rows, n_cols + , eigenval_begin, eigenvec_begin) ; + for (unsigned int lane = 1; lane < macro_size; ++lane) + { + std::advance (mass_cbegin, nm) ; + std::advance (deriv_cbegin, nm) ; + std::advance (eigenvec_begin, nm) ; + std::advance (eigenval_begin, n_rows) ; + spectral_assembly (mass_cbegin, deriv_cbegin, n_rows, n_cols + , eigenval_begin, eigenvec_begin) ; + } + + eigenvalues[dir].resize (n_rows) ; + 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)) ; +} + + + +template +inline +void +TensorProductMatrixSymmetricSum,size> +::reinit (const std::array >,dim> &mass_matrix, + const std::array >,dim> &derivative_matrix) +{ + reinit_impl (mass_matrix, derivative_matrix) ; +} + + + +template +inline +void +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())); + TensorProductMatrixSymmetricSumBase,size>::vmult (dst.begin(), src.begin()); +} + + + +template +inline +void +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())); + TensorProductMatrixSymmetricSumBase,size>::apply_inverse (dst.begin(), src.begin()); }