#include <deal.II/lac/tensor_product_matrix.h>

Public Types
using	value_type = Number

Public Member Functions
	TensorProductMatrixSymmetricSum ()=default

template<typename T >
	TensorProductMatrixSymmetricSum (const T &mass_matrix, const T &derivative_matrix)

template<typename T >
void	reinit (const T &mass_matrix, const T &derivative_matrix)

unsigned int	m () const

unsigned int	n () const

void	vmult (const ArrayView< Number > &dst, const ArrayView< const Number > &src) const

void	vmult (const ArrayView< Number > &dst, const ArrayView< const Number > &src, AlignedVector< Number > &tmp) const

void	apply_inverse (const ArrayView< Number > &dst, const ArrayView< const Number > &src) const

std::size_t	memory_consumption () const

Static Public Attributes
static constexpr int	n_rows_1d_static = n_rows_1d

Protected Attributes
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

Private Attributes
AlignedVector< Number >	tmp_array

Threads::Mutex	mutex

Detailed Description

template<int dim, typename Number, int n_rows_1d = -1>
class TensorProductMatrixSymmetricSum< dim, Number, n_rows_1d >

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 1964's work by Lynch et al. [126],

\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: 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 n_rows_1d 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.; This class can work with scalar types (float, double) and VectorizedArray types.

Template Parameters

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.
n_rows_1d	Compile-time number of rows of 1d matrices (only valid if the number of rows and columns coincide for each dimension). By default at -1, which means that the number of rows is determined at run-time by means of the matrices passed to the reinit() function.

Definition at line 113 of file tensor_product_matrix.h.