Home Installation C++ docs Python docs
Functions
basix::precompute Namespace Reference

Functions

std::vector< std::size_t > prepare_permutation (const std::vector< std::size_t > &perm)
 
template<typename E >
void apply_permutation (const std::vector< std::size_t > &perm, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
 
template<typename E >
void apply_permutation_to_transpose (const std::vector< std::size_t > &perm, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
 
std::tuple< std::vector< std::size_t >, std::vector< double >, xt::xtensor< double, 2 > > prepare_matrix (const xt::xtensor< double, 2 > &matrix)
 
template<typename T , typename E >
void apply_matrix (const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
 
template<typename T , typename E >
void apply_matrix_to_transpose (const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &matrix, const xtl::span< E > &data, std::size_t offset=0, std::size_t block_size=1)
 

Detailed Description

Matrix and permutation precomputation

These functions generate precomputed version of matrices to allow application without temporary memory assignment later

Function Documentation

◆ apply_matrix()

template<typename T , typename E >
void basix::precompute::apply_matrix ( const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &  matrix,
const xtl::span< E > &  data,
std::size_t  offset = 0,
std::size_t  block_size = 1 
)

Apply a (precomputed) matrix

This uses the representation returned by prepare_matrix() to apply a matrix without needing any temporary memory.

In pseudo code, this function does the following:

perm, diag, mat = matrix
apply_permutation(perm, data)
FOR index IN RANGE(dim):
data[index] *= diag[index]
FOR j IN RANGE(dim):
data[index] *= mat[index, j] * data[j]

If block_size is set, this will apply the permutation to every block. The offset is set, this will start applying the permutation at the offsetth block.

Example

As an example, consider the matrix \(A = \) [[-1, 0, 1], [1, 1, 0], [2, 0, 2]]. In the documentation of prepare_matrix(), we saw that the precomputed representation of this matrix is the identity permutation,

\begin{align*} D &= \begin{bmatrix}-1\\1\\4\end{bmatrix},\\ \quad M &= \begin{bmatrix} 0&0&1\\ -1&0&1\\ -2&0&0 \end{bmatrix}. \end{align*}

In this example, we look at how this representation can be used to apply this matrix to the vector \(v = \) [3, -1, 2].

No permutation is necessary, so first, we multiply \(v_0\) by \(D_0=-1\). After this, \(v\) is [-3, -1, 2].

Next, we add \(M_{0,i}v_i\) to \(v_0\) for all \(i\): in this case, we add \(0\times-3 + 0\times-1 + 1\times2 = 2\). After this, \(v\) is [-1, -1, 2].

Next, we multiply \(v_1\) by \(D_1=1\). After this, \(v\) is [-1, -1, 2].

Next, we add \(M_{1,i}v_i\) to \(v_1\) for all \(i\): in this case, we add \(-1\times-1 + 0\times-1 + 1\times2 = 3\). After this, \(v\) is [-1, 2, 2].

Next, we multiply \(v_2\) by \(D_2=4\). After this, \(v\) is [-1, 2, 8].

Next, we add \(M_{2,i}v_i\) to \(v_2\) for all \(i\): in this case, we add \(-2\times-1 + 0\times2 + 0\times8 = 2\). After this, \(v\) is [-1, 2, 10]. This final value of \(v\) is what the result of \(Av\)

Parameters
[in]matrixA matrix in precomputed form (as returned by prepare_matrix())
[in,out]dataThe data to apply the permutation to
[in]offsetThe position in the data to start applying the permutation
[in]block_sizeThe block size of the data

◆ apply_matrix_to_transpose()

template<typename T , typename E >
void basix::precompute::apply_matrix_to_transpose ( const std::tuple< std::vector< std::size_t >, std::vector< T >, xt::xtensor< T, 2 >> &  matrix,
const xtl::span< E > &  data,
std::size_t  offset = 0,
std::size_t  block_size = 1 
)

Apply a (precomputed) matrix to some transposed data.

See apply_matrix().

◆ apply_permutation()

template<typename E >
void basix::precompute::apply_permutation ( const std::vector< std::size_t > &  perm,
const xtl::span< E > &  data,
std::size_t  offset = 0,
std::size_t  block_size = 1 
)

Apply a (precomputed) permutation

This uses the representation returned by prepare_permutation() to apply a permutation without needing any temporary memory.

In pseudo code, this function does the following:

FOR index, entry IN perm:
SWAP(data[index], data[entry])

If block_size is set, this will apply the permutation to every block. The offset is set, this will start applying the permutation at the offsetth block.

Example

As an example, consider the permutation P = [1, 4, 0, 5, 2, 3]. In the documentation of prepare_permutation(), we saw that the precomputed representation of this permutation is P2 = [1, 4, 4, 5, 4, 5]. In this example, we look at how this representation can be used to apply this permutation to the array A = [a, b, c, d, e, f].

P2[0] is 1, so we swap A[0] and A[1]. After this, A = [b, a, c, d, e, f].

P2[1] is 4, so we swap A[1] and A[4]. After this, A = [b, e, c, d, a, f].

P2[2] is 4, so we swap A[2] and A[4]. After this, A = [b, e, a, d, c, f].

P2[3] is 5, so we swap A[3] and A[5]. After this, A = [b, e, a, f, c, d].

P2[4] is 4, so we swap A[4] and A[4]. This changes nothing.

P2[5] is 5, so we swap A[5] and A[5]. This changes nothing.

Therefore the result of applying this permutation is [b, e, a, f, c, d] (which is what we get if we apply the permutation directly).

Parameters
[in]permA permutation in precomputed form (as returned by prepare_permutation())
[in,out]dataThe data to apply the permutation to
[in]offsetThe position in the data to start applying the permutation
[in]block_sizeThe block size of the data

◆ apply_permutation_to_transpose()

template<typename E >
void basix::precompute::apply_permutation_to_transpose ( const std::vector< std::size_t > &  perm,
const xtl::span< E > &  data,
std::size_t  offset = 0,
std::size_t  block_size = 1 
)

Apply a (precomputed) permutation to some transposed data

see apply_permutation().

◆ prepare_matrix()

std::tuple< std::vector< std::size_t >, std::vector< double >, xt::xtensor< double, 2 > > basix::precompute::prepare_matrix ( const xt::xtensor< double, 2 > &  matrix)

Prepare a matrix

This computes a representation of the matrix that allows the matrix to be applied without any temporary memory assignment.

This function will first permute the matrix's columns so that the top left \(n\times n\) blocks are invertible (for all \(n\)). Let \(A\) be the input matrix after the permutation is applied. The output vector \(D\) and matrix \(M\) are then given by:

\begin{align*} D_i &= \begin{cases} A_{i, i} & i = 0\\ A_{i, i} - A_{i,:i}A_{:i,:i}^{-1}A_{:i,i} & i \not= 0 \end{cases},\\ M_{i,j} &= \begin{cases} A_{i,:i}A_{:i,:i}^{-1}e_j & j < i\\ 0 & j = i\\ A_{i, i} - A_{i,:i}A_{:i,:i}^{-1}A_{:i,j} & j > i = 0 \end{cases}, \end{align*}

where \(e_j\) is the \(j\)th coordinate vector, we index all the matrices and vector starting at 0, and we use numpy-slicing-stying notation in the subscripts: for example, \(A_{:i,j}\) represents the first \(i\) entries in the \(j\)th column of \(A\)

This function returns the permutation (precomputed as in prepare_permutation()), the vector \(D\), and the matrix \(M\) as a tuple.

Example

As an example, consider the matrix \(A = \) [[-1, 0, 1], [1, 1, 0], [2, 0, 2]]. For this matrix, no permutation is needed, so the first item in the output will represent the identity permutation. We now compute the output vector \(D\) and matrix \(M\).

First, we set \(D_0 = A_{0,0}=-1\), set the diagonal of \(M\) to be 0 and set \(M_{0, 1:} = A_{0, 1:}=\begin{bmatrix}0&1\end{bmatrix}\). The output so far is

\begin{align*} D &= \begin{bmatrix}-1\\?\\?\end{bmatrix},\\ \quad M &= \begin{bmatrix} 0&0&1\\ ?&0&?\\ ?&?&0 \end{bmatrix}. \end{align*}

Next, we set:

\begin{align*} D_1 &= A_{1,1} - A_{1, :1}A_{:1,:1}^{-1}A_{:1, 1}\\ &= 1 - \begin{bmatrix}-1\end{bmatrix}\cdot\begin{bmatrix}0\end{bmatrix}\\ &= 1,\\ M_{2,0} &= A_{1, :1}A_{:1,:1}^{-1}e_0\\ &= \begin{bmatrix}1\end{bmatrix}\begin{bmatrix}-1\end{bmatrix}^{-1} \begin{bmatrix}1\end{bmatrix}\\ &= \begin{bmatrix}-1\end{bmatrix} M_{2,3} &= A_{1,2}-A_{1, :1}A_{:1,:1}^{-1}A_{:1, 1}\\ &= 0-\begin{bmatrix}1\end{bmatrix}\begin{bmatrix}-1\end{bmatrix}^{-1} \begin{bmatrix}1\end{bmatrix},\\ &= 1. \end{align*}

The output so far is

\begin{align*} D &= \begin{bmatrix}-1\\1\\?\end{bmatrix},\\ \quad M &= \begin{bmatrix} 0&0&1\\ -1&0&1\\ ?&?&0 \end{bmatrix}. \end{align*}

Next, we set:

\begin{align*} D_2 &= A_{2,2} - A_{2, :2}A_{:2,:2}^{-1}A_{:2, 2}\\ &= 2 - \begin{bmatrix}2&0\end{bmatrix} \begin{bmatrix}-1&0\\1&1\end{bmatrix}^{-1} \begin{bmatrix}1\\0\end{bmatrix}\\ &= 4,\\ M_{2,0} &= A_{2, :2}A_{:2,:2}^{-1}e_0\\ &= -2.\\ M_{2,1} &= A_{2, :2}A_{:2,:2}^{-1}e_1\\ &= 0.\\ \end{align*}

The output is

\begin{align*} D &= \begin{bmatrix}-1\\1\\4\end{bmatrix},\\ \quad M &= \begin{bmatrix} 0&0&1\\ -1&0&1\\ -2&0&0 \end{bmatrix}. \end{align*}

For an example of how the permutation in this form is applied, see apply_matrix().

Parameters
[in]matrixA matrix
Returns
The precomputed representation of the matrix

◆ prepare_permutation()

std::vector< std::size_t > basix::precompute::prepare_permutation ( const std::vector< std::size_t > &  perm)

Prepare a permutation

This computes a representation of the permutation that allows the permutation to be applied without any temporary memory assignment.

In pseudo code, this function does the following:

FOR index, entry IN perm:
new_index = entry
WHILE new_index < index:
new_index = perm[new_index]
OUT[index] = new_index

Example

As an example, consider the permutation P = [1, 4, 0, 5, 2, 3].

First, we look at the 0th entry. P[0] is 1. This is greater than 0, so the 0th entry of the output is 1.

Next, we look at the 1st entry. P[1] is 4. This is greater than 1, so the 1st entry of the output is 4.

Next, we look at the 2nd entry. P[2] is 0. This is less than 2, so we look at P[0].P[0]is 1. This is less than 2, so we look atP[1].P[1]` is 4. This is greater than 2, so the 2nd entry of the output is 4.

Next, we look at the 3rd entry. P[3] is 5. This is greater than 3, so the 3rd entry of the output is 5.

Next, we look at the 4th entry. P[4] is 2. This is less than 4, so we look at P[2]. P[2] is 0. This is less than 4, so we look at P[0]. P[0] is 1. This is less than 4, so we look at P[1]. P[1] is 4. This is greater than (or equal to) 4, so the 4th entry of the output is 4.

Next, we look at the 5th entry. P[5] is 3. This is less than 5, so we look at P[3]. P[3] is 5. This is greater than (or equal to) 5, so the 5th entry of the output is 5.

Hence, the output of this function in this case is [1, 4, 4, 5, 4, 5].

For an example of how the permutation in this form is applied, see apply_permutation().

Parameters
[in]permA permutation
Returns
The precomputed representation of the permutation