![]() |
Reference documentation for deal.II version 9.1.1
|
#include <deal.II/lac/tensor_product_matrix.h>
Public Member Functions | |
TensorProductMatrixSymmetricSum ()=default | |
TensorProductMatrixSymmetricSum (const std::array< Table< 2, Number >, dim > &mass_matrix, const std::array< Table< 2, Number >, dim > &derivative_matrix) | |
TensorProductMatrixSymmetricSum (const std::array< FullMatrix< Number >, dim > &mass_matrix, const std::array< FullMatrix< Number >, dim > &derivative_matrix) | |
TensorProductMatrixSymmetricSum (const Table< 2, Number > &mass_matrix, const Table< 2, Number > &derivative_matrix) | |
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) |
void | reinit (const Table< 2, Number > &mass_matrix, const Table< 2, Number > &derivative_matrix) |
![]() | |
unsigned int | m () const |
unsigned int | n () const |
void | vmult (const ArrayView< Number > &dst, const ArrayView< const Number > &src) const |
void | apply_inverse (const ArrayView< Number > &dst, const ArrayView< const Number > &src) const |
Private Member Functions | |
template<typename MatrixArray > | |
void | reinit_impl (MatrixArray &&mass_matrix, MatrixArray &&derivative_matrix) |
Additional Inherited Members | |
![]() | |
TensorProductMatrixSymmetricSumBase ()=default | |
![]() | |
std::array< Table< 2, Number >, dim > | mass_matrix |
std::array< Table< 2, Number >, dim > | derivative_matrix |
std::array< AlignedVector< Number >, dim > | eigenvalues |
std::array< Table< 2, Number >, dim > | eigenvectors |
This is a special matrix class defined as the tensor product (or Kronecker product) of 1D matrices of the type
\begin{align*} L &= A_1 \otimes M_0 + M_1 \otimes A_0 \end{align*}
in 2D and
\begin{align*} L &= A_2 \otimes M_1 \otimes M_0 + M_2 \otimes A_1 \otimes M_0 + M_2 \otimes M_1 \otimes A_0 \end{align*}
in 3D. The typical application setting is a discretization of the Laplacian \(L\) on a Cartesian (axis-aligned) geometry, where it can be exactly represented by the Kronecker or tensor product of a 1D mass matrix \(M\) and a 1D Laplace matrix \(A\) in each tensor direction (due to symmetry \(M\) and \(A\) are the same in each dimension). The dimension of the resulting class is the product of the one-dimensional matrices.
This class implements two basic operations, namely the usual multiplication by a vector and the inverse. For both operations, fast tensorial techniques can be applied that implement the operator evaluation in \(\text{size}(M)^{d+1}\) arithmetic operations, considerably less than \(\text{size}(M)^{2d}\) for the naive forward transformation and \(\text{size}(M)^{3d}\) for setting up the inverse of \(L\).
Interestingly, the exact inverse of the matrix \(L\) can be found through tensor products due to an article by R. E. Lynch, J. R. Rice, D. H. Thomas, Direct solution of partial difference equations by tensor product methods, Numerische Mathematik 6, 185-199 from 1964,
\begin{align*} L^{-1} &= S_1 \otimes S_0 (\Lambda_1 \otimes I + I \otimes \Lambda_0)^{-1} S_1^\mathrm T \otimes S_0^\mathrm T, \end{align*}
where \(S_d\) is the matrix of eigenvectors to the generalized eigenvalue problem in the given tensor direction \(d\):
\begin{align*} A_d s &= \lambda M_d s, d = 0, \quad \ldots,\mathrm{dim}, \end{align*}
and \(\Lambda_d\) is the diagonal matrix representing the generalized eigenvalues \(\lambda\). Note that the vectors \(s\) are such that they simultaneously diagonalize \(A_d\) and \(M_d\), i.e. \(S_d^{\mathrm T} A_d S_d = \Lambda_d\) and \(S_d^{\mathrm T} M_d S_d = I\). This method of matrix inversion is called fast diagonalization method.
This class requires LAPACK support.
Note that this class allows for two modes of usage. The first is a use case with run time constants for the matrix dimensions that is achieved by setting the optional template parameter for the size to -1. The second mode of usage that is faster allows to set the template parameter as a compile time constant, giving significantly faster code in particular for small sizes of the matrix.
dim | Dimension of the problem. Currently, 1D, 2D, and 3D codes are implemented. |
Number | Arithmetic type of the underlying array elements. Note that the underlying LAPACK implementation supports only float and double numbers, so only these two types are currently supported by the generic class. Nevertheless, a template specialization for the vectorized types VectorizedArray<float> and VectorizedArray<double> exists. This is necessary to perform LAPACK calculations for each vectorization lane, i.e. for the supported float and double numbers. |
size | Compile-time array lengths. By default at -1, which means that the run-time information stored in the matrices passed to the reinit() function is used. |
Definition at line 225 of file tensor_product_matrix.h.
|
default |
Default constructor.
TensorProductMatrixSymmetricSum< dim, Number, size >::TensorProductMatrixSymmetricSum | ( | const std::array< Table< 2, Number >, dim > & | mass_matrix, |
const std::array< Table< 2, Number >, dim > & | derivative_matrix | ||
) |
Constructor that is equivalent to the empty constructor and immediately calling reinit(const std::array<Table<2,Number>, dim>&,const std::array<Table<2,Number>, dim>&).
TensorProductMatrixSymmetricSum< dim, Number, size >::TensorProductMatrixSymmetricSum | ( | const std::array< FullMatrix< Number >, dim > & | mass_matrix, |
const std::array< FullMatrix< Number >, dim > & | derivative_matrix | ||
) |
Constructor that is equivalent to the empty constructor and immediately calling reinit(const std::array<FullMatrix<Number>,dim>&,const std::array<FullMatrix<Number>,dim>&).
TensorProductMatrixSymmetricSum< dim, Number, size >::TensorProductMatrixSymmetricSum | ( | const Table< 2, Number > & | mass_matrix, |
const Table< 2, Number > & | derivative_matrix | ||
) |
Constructor that is equivalent to the empty constructor and immediately calling reinit(const Table<2,Number>&,const Table<2,Number>&).
void TensorProductMatrixSymmetricSum< dim, Number, size >::reinit | ( | const std::array< Table< 2, Number >, dim > & | mass_matrix, |
const std::array< Table< 2, Number >, dim > & | derivative_matrix | ||
) |
Initializes the tensor product matrix by copying the arrays of 1D mass matrices mass_matrix
and 1D derivative matrices derivative_matrix
into its base class counterparts, respectively, and by assembling the regarding generalized eigenvalues and eigenvectors in TensorProductMatrixSymmetricSumBase::eigenvalues and TensorProductMatrixSymmetricSumBase::eigenvectors, respectively. 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 definite.
void TensorProductMatrixSymmetricSum< dim, Number, size >::reinit | ( | const std::array< FullMatrix< Number >, dim > & | mass_matrix, |
const std::array< FullMatrix< Number >, dim > & | derivative_matrix | ||
) |
This function is equivalent to the previous reinit() except that the 1D matrices in mass_matrix
and derivative_matrix
are passed in terms of a FullMatrix, respectively.
void TensorProductMatrixSymmetricSum< dim, Number, size >::reinit | ( | const Table< 2, Number > & | mass_matrix, |
const Table< 2, Number > & | derivative_matrix | ||
) |
This function is equivalent to the first reinit() except that we consider the same 1D mass matrix mass_matrix
and the same 1D derivative matrix derivative_matrix
for each tensor direction.
|
private |
A generic implementation of all reinit() functions based on perfect forwarding, that allows to pass lvalue as well as rvalue arguments.
MatrixArray | Has to be convertible to the underlying type of TensorProductMatrixSymmetricSumBase::mass_matrix and TensorProductMatrixSymmetricSumBase::derivative_matrix. |