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)

## Matrix and permutation precomputation

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

## ◆ 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] matrix A matrix in precomputed form (as returned by prepare_matrix()) [in,out] data The data to apply the permutation to [in] offset The position in the data to start applying the permutation [in] block_size The 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] perm A permutation in precomputed form (as returned by prepare_permutation()) [in,out] data The data to apply the permutation to [in] offset The position in the data to start applying the permutation [in] block_size The 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] matrix A 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] perm A permutation
Returns
The precomputed representation of the permutation