Fixed-size objects

tmatrix and tvector: the dimension is read from their indexing policy at compile time, and they resolve the compile-time solver with the internal residency.

Factorize, then substitute

One factorization, two right-hand sides.

\[\begin{split} A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right), \quad b_1 = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad b_2 = \begin{pmatrix} 7 \\ 7 \\ 8 \\ 10 \end{pmatrix}, \quad x_1 = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix}, \quad x_2 = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
const tfel::math::tvector<4, double> b1{14, 14, 24, 33}, b2{7, 7, 8, 10};
tfel::math::tvector<4, double> x1, x2;
tfel::math::fsarray<4, int> piv;

// the dimension, the residencies and the strides are deduced from the objects
if (!tdls::factorize(A, piv)) return 1;

tdls::substitute(A, piv, b1, x1);
tdls::substitute(A, piv, b2, x2);

Factorize, then substitute in place

One buffer for the right-hand side and the solution.

\[\begin{split} A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tvector<4, double> y{14, 14, 24, 33};
tfel::math::fsarray<4, int> piv;

if (!tdls::factorize(A, piv)) return 1;

// y holds b on entry and x on exit
tdls::substitute_inplace(A, piv, y);

Solve in one call

Factorization and substitution in one call, two buffers.

\[\begin{split} A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
const tfel::math::tvector<4, double> b{14, 14, 24, 33};
tfel::math::tvector<4, double> x;
tfel::math::fsarray<4, int> piv;

const bool ok = tdls::solve(A, piv, b, x);

Solve in place

Factorization and substitution in one call, one buffer.

\[\begin{split} A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tvector<4, double> y{14, 14, 24, 33};
tfel::math::fsarray<4, int> piv;

// y holds b on entry and x on exit
const bool ok = tdls::solve_inplace(A, piv, y);

One canonical column

The right-hand side is a canonical vector, generated on the fly.

\[\begin{split} A x = e_2, \quad A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right), \quad e_2 = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tvector<4, double> x;
tfel::math::fsarray<4, int> piv;

if (!tdls::factorize(A, piv)) return 1;

// x is column 2 of A^-1
tdls::substitute_canonical(A, piv, 2, x);

With a configuration value

A configuration value of the matrix scalar type, as the first template argument. Every entry point has this form. Tiles of 2 here.

\[\begin{split} A = \left(\begin{array}{cc|cc} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ \hline 0 & 1 & 6 & 1 \\ 2 & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tvector<4, double> y{14, 14, 24, 33};
tfel::math::fsarray<4, int> piv;

// the configuration value of the matrix scalar type, as first template argument
constexpr auto config =
    tdls::TiledLUppConfig<double>{.tile_size = 2, .schedule = tdls::Schedule::LeftLooking};
const bool ok = tdls::solve_inplace<config>(A, piv, y);

Counting overloads

factorize, solve and solve_inplace take a trailing int&: the number of columns whose best in-tile pivot fell below the threshold. The three candidates of column 0 inside the tile, in red, are tiny; the search below the tile takes row 3, in blue.

\[\begin{split} A = \left(\begin{array}{ccc|c} \textcolor{red}{10^{-12}} & 1 & 0 & 2 \\ \textcolor{red}{2 \cdot 10^{-12}} & 5 & 1 & 0 \\ \textcolor{red}{3 \cdot 10^{-12}} & 1 & 6 & 1 \\ \hline \textcolor{blue}{2} & 0 & 1 & 7 \end{array}\right), \quad b = \begin{pmatrix} 10 \\ 13 \\ 24 \\ 31 \end{pmatrix}, \quad x = \begin{pmatrix} 0 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

tfel::math::tmatrix<4, 4, double> A{1e-12, 1, 0, 2, 2e-12, 5, 1, 0, 3e-12, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tvector<4, double> y{10, 13, 24, 31};
tfel::math::fsarray<4, int> piv;

// the trailing int&: the columns whose best in-tile pivot fell below the threshold
int oot_count;
const bool ok = tdls::solve_inplace(A, piv, y, oot_count);

const bool searched_below_the_tile = oot_count == 1 && piv[0] == 3;

Pivot forms

The pivot is any dense int object of the dimension, or a raw int array or pointer. The four factorizations are bitwise identical.

\[\begin{split} A = \left(\begin{array}{ccc|c} 4 & 1 & 0 & 2 \\ 1 & 5 & 1 & 0 \\ 0 & 1 & 6 & 1 \\ \hline 2 & 0 & 1 & 7 \end{array}\right) \end{split}\]

tfel::math::tmatrix<4, 4, double> A1{4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
tfel::math::tmatrix<4, 4, double> A2 = A1, A3 = A1, A4 = A1;

// the pivot is any dense int object of extent 4, or a raw int array or pointer
tfel::math::fsarray<4, int> piv_fsarray;
tfel::math::tvector<4, int> piv_tvector;
int piv_array[4];
int storage[4];
int* piv_pointer = storage;

const bool ok = tdls::factorize(A1, piv_fsarray) && tdls::factorize(A2, piv_tvector) &&
                tdls::factorize(A3, piv_array) && tdls::factorize(A4, piv_pointer);