Configuration
The knobs of TiledLUppConfig, one snippet each. Their roles are
listed on the family page.
Defaults
The default value, and the two solvers named on it. Tiles of 3 cut this 4 x 4 matrix into a full tile and a trailing tile of 1.
\[\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}\]
// every knob at its default: tiles of 3, right-looking, default thresholds
constexpr tdls::TiledLUppConfig<double> config{};
using Static = tdls::TiledLUppSolverStatic<double, 4, config>;
using Dynamic = tdls::TiledLUppSolverDynamic<double, config>;
// the configuration argument of both solvers defaults to that value
static_assert(std::is_same_v<Static, tdls::TiledLUppSolverStatic<double, 4>>);
static_assert(std::is_same_v<Dynamic, tdls::TiledLUppSolverDynamic<double>>);
int piv[4];
// y holds b on entry and x on exit
const bool ok_static = Static::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);
const bool ok_dynamic = Dynamic::solve_inplace(4, A2, 1, piv, 1, y2, 1);
Tile size
The same 6 x 6 matrix under two tile sizes: one full tile of 4 with a trailing tile of 2, or a 2 x 2 grid of full tiles of 3.
\[\begin{split}
\left(\begin{array}{cccc|cc}
8 & 1 & 0 & 0 & 2 & 0 \\
1 & 8 & 1 & 0 & 0 & 2 \\
0 & 1 & 8 & 1 & 0 & 0 \\
0 & 0 & 1 & 8 & 1 & 0 \\ \hline
2 & 0 & 0 & 1 & 8 & 1 \\
0 & 2 & 0 & 0 & 1 & 8
\end{array}\right)
\qquad
\left(\begin{array}{ccc|ccc}
8 & 1 & 0 & 0 & 2 & 0 \\
1 & 8 & 1 & 0 & 0 & 2 \\
0 & 1 & 8 & 1 & 0 & 0 \\ \hline
0 & 0 & 1 & 8 & 1 & 0 \\
2 & 0 & 0 & 1 & 8 & 1 \\
0 & 2 & 0 & 0 & 1 & 8
\end{array}\right)
\end{split}\]
\[\begin{split}
b = \begin{pmatrix} 20 \\ 32 \\ 30 \\ 40 \\ 52 \\ 57 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{pmatrix}
\end{split}\]
constexpr auto tiles_of_4 = tdls::TiledLUppConfig<double>{.tile_size = 4};
constexpr auto tiles_of_3 = tdls::TiledLUppConfig<double>{.tile_size = 3};
// 6 x 6: one full tile of 4 and a trailing tile of 2, or a 2 x 2 grid of full tiles
using Tiles4 = tdls::TiledLUppSolverStatic<double, 6, tiles_of_4>;
using Tiles3 = tdls::TiledLUppSolverStatic<double, 6, tiles_of_3>;
static_assert(Tiles4::full_tiles == 1 && Tiles4::last_tile_tail == 2);
static_assert(Tiles3::full_tiles == 2 && Tiles3::last_tile_tail == 0);
int piv[6];
const bool ok4 = Tiles4::solve_inplace<true, true, true>(A4, 1, piv, 1, y4, 1);
const bool ok3 = Tiles3::solve_inplace<true, true, true>(A3, 1, piv, 1, y3, 1);
A single tile
A tile size equal to the dimension: the whole matrix is one tile.
\[\begin{split}
A = \begin{pmatrix}
4 & 1 & 0 & 2 \\
1 & 5 & 1 & 0 \\
0 & 1 & 6 & 1 \\
2 & 0 & 1 & 7
\end{pmatrix}, \quad
b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix}
\end{split}\]
// one tile holds the whole matrix: no trailing update at all
constexpr auto one_tile = tdls::TiledLUppConfig<double>{.tile_size = 4};
using Solver = tdls::TiledLUppSolverStatic<double, 4, one_tile>;
static_assert(Solver::full_tiles == 1 && Solver::last_tile_tail == 0);
int piv[4];
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);
A tile larger than the dimension
A tile size above the dimension: the grid is a single partial tile.
\[\begin{split}
A = \begin{pmatrix}
3 & 1 & 0 \\
1 & 4 & 1 \\
0 & 1 & 5
\end{pmatrix}, \quad
b = \begin{pmatrix} 5 \\ 12 \\ 17 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix}
\end{split}\]
// the grid is a single partial tile: 3 of its 8 x 8 slots are used, the others never touched
constexpr auto wide = tdls::TiledLUppConfig<double>{.tile_size = 8};
using Solver = tdls::TiledLUppSolverStatic<double, 3, wide>;
static_assert(Solver::full_tiles == 0 && Solver::last_tile_tail == 3);
int piv[3];
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);
Scalar tiles
Tiles of one element: every step handles one pivot.
\[\begin{split}
A = \left(\begin{array}{c|c|c|c}
4 & 1 & 0 & 2 \\ \hline
1 & 5 & 1 & 0 \\ \hline
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}\]
// tiles of one element: an untiled scalar elimination
constexpr auto scalar = tdls::TiledLUppConfig<double>{.tile_size = 1};
using Solver = tdls::TiledLUppSolverStatic<double, 4, scalar>;
static_assert(Solver::full_tiles == 4 && Solver::last_tile_tail == 0);
int piv[4];
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);
Schedules
The two elimination schedules on the same system.
\[\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}\]
constexpr auto right =
tdls::TiledLUppConfig<double>{.tile_size = 2, .schedule = tdls::Schedule::RightLooking};
constexpr auto left =
tdls::TiledLUppConfig<double>{.tile_size = 2, .schedule = tdls::Schedule::LeftLooking};
using RightLooking = tdls::TiledLUppSolverStatic<double, 4, right>;
using LeftLooking = tdls::TiledLUppSolverStatic<double, 4, left>;
int piv_right[4], piv_left[4];
// the same factors from two update orders
const bool ok_right =
RightLooking::solve_inplace<true, true, true>(A_right, 1, piv_right, 1, y_right, 1);
const bool ok_left =
LeftLooking::solve_inplace<true, true, true>(A_left, 1, piv_left, 1, y_left, 1);
Pivot thresholds
The last pivot is tiny. The default thresholds accept it, raised ones declare the matrix singular.
\[\begin{split}
A = \left(\begin{array}{cc|cc}
4 & 1 & 0 & 0 \\
1 & 5 & 0 & 0 \\ \hline
0 & 0 & 6 & 0 \\
0 & 0 & 0 & \textcolor{red}{10^{-8}}
\end{array}\right), \quad
b = \begin{pmatrix} 6 \\ 11 \\ 18 \\ 0 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 0 \end{pmatrix}
\end{split}\]
constexpr auto lenient = tdls::TiledLUppConfig<double>{.tile_size = 2};
constexpr auto strict = tdls::TiledLUppConfig<double>{
.tile_size = 2, .oot_threshold = 1e-6, .singular_floor = 1e-6};
using Lenient = tdls::TiledLUppSolverStatic<double, 4, lenient>;
using Strict = tdls::TiledLUppSolverStatic<double, 4, strict>;
int piv[4];
// the last pivot, 1e-8, reaches the default threshold: accepted inside its tile
const bool ok_lenient = Lenient::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);
// below 1e-6 the search finds no row left, and the floor declares the matrix singular
const bool ok_strict = Strict::solve_inplace<true, true, true>(A2, 1, piv, 1, y2, 1);
Full out-of-tile scan
Column 0 has no acceptable pivot inside its tile. Two rows below can provide one: the first acceptable, in red, or the largest, in blue.
\[\begin{split}
A = \left(\begin{array}{cc|cc}
10^{-12} & 1 & 0 & 2 \\
2 \cdot 10^{-12} & 5 & 1 & 0 \\ \hline
\textcolor{red}{10^{-6}} & 1 & 6 & 1 \\
\textcolor{blue}{3} & 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}\]
constexpr auto first = tdls::TiledLUppConfig<double>{.tile_size = 2};
constexpr auto best =
tdls::TiledLUppConfig<double>{.tile_size = 2, .oot_first_acceptable = false};
using First = tdls::TiledLUppSolverStatic<double, 4, first>;
using Best = tdls::TiledLUppSolverStatic<double, 4, best>;
int piv_first[4], piv_best[4];
// both in-tile candidates of column 0 are below 1e-10: the search goes below the tile
const bool ok_first = First::factorize<true, true>(A_first, 1, piv_first, 1);
const bool ok_best = Best::factorize<true, true>(A_best, 1, piv_best, 1);
// first acceptable: row 2, whose 1e-6 reaches the threshold; full scan: row 3, the largest
const bool rows = piv_first[0] == 2 && piv_best[0] == 3;
No forced unrolling
The same solve with and without the unroll pragmas: bitwise identical factors and solution.
\[\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}\]
constexpr auto unrolled = tdls::TiledLUppConfig<double>{.tile_size = 2};
constexpr auto rolled = tdls::TiledLUppConfig<double>{.tile_size = 2, .unroll_inner = false};
using Unrolled = tdls::TiledLUppSolverStatic<double, 4, unrolled>;
using Rolled = tdls::TiledLUppSolverStatic<double, 4, rolled>;
int piv[4];
// no unroll pragma: same values, faster builds
const bool ok_unrolled =
Unrolled::solve_inplace<true, true, true>(A_unrolled, 1, piv, 1, y_unrolled, 1);
const bool ok_rolled =
Rolled::solve_inplace<true, true, true>(A_rolled, 1, piv, 1, y_rolled, 1);
Column-major layout
The matrix stored column by column, as one system of a batch of 3: the layout knob remaps the flat index, the stride does the rest.
\[\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)
\ \text{stored as}\
\begin{pmatrix}
4 & 1 & 0 & 2 \\
1 & 5 & 1 & 0 \\
0 & 1 & 6 & 1 \\
2 & 0 & 1 & 7
\end{pmatrix}^{\mathsf T}, \quad
b = \begin{pmatrix} 14 \\ 14 \\ 24 \\ 33 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix}
\end{split}\]
constexpr auto colmajor =
tdls::TiledLUppConfig<double>{.tile_size = 2, .layout = tdls::MatrixLayout::ColMajor};
using Solver = tdls::TiledLUppSolverStatic<double, 4, colmajor>;
// system 1 of an SoA batch of 3: element k of the matrix sits at batch[k * 3 + 1],
// with k = c * 4 + r under the column-major layout
double* A = batch + 1;
const int stride = 3;
int piv[4];
const bool ok = Solver::solve_inplace<true, true, false>(A, stride, piv, 1, y, 1);
Scalar types
The same system in float and in long double.
\[\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}\]
// float: the acceptable-pivot threshold defaults to 1e-4
using Single = tdls::TiledLUppSolverStatic<float, 4>;
static_assert(Single::oot_threshold == 1e-4f);
// long double: the thresholds are written as double literals and stored exactly
constexpr auto extended =
tdls::TiledLUppConfig<long double>{.oot_threshold = 1e-12, .singular_floor = 1e-300};
using Extended = tdls::TiledLUppSolverStatic<long double, 4, extended>;
static_assert(Extended::oot_threshold == static_cast<long double>(1e-12));
int piv[4];
const bool ok_single =
Single::solve_inplace<true, true, true>(A_single, 1, piv, 1, y_single, 1);
const bool ok_extended =
Extended::solve_inplace<true, true, true>(A_extended, 1, piv, 1, y_extended, 1);
Solver constants
The tile grid and the knobs, read back from the solver type.
\[\begin{split}
A = \left(\begin{array}{cccc|cc}
8 & 1 & 0 & 0 & 2 & 0 \\
1 & 8 & 1 & 0 & 0 & 2 \\
0 & 1 & 8 & 1 & 0 & 0 \\
0 & 0 & 1 & 8 & 1 & 0 \\ \hline
2 & 0 & 0 & 1 & 8 & 1 \\
0 & 2 & 0 & 0 & 1 & 8
\end{array}\right), \quad
b = \begin{pmatrix} 20 \\ 32 \\ 30 \\ 40 \\ 52 \\ 57 \end{pmatrix}, \quad
x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{pmatrix}
\end{split}\]
constexpr auto config =
tdls::TiledLUppConfig<double>{.tile_size = 4, .schedule = tdls::Schedule::LeftLooking};
using Solver = tdls::TiledLUppSolverStatic<double, 6, config>;
// the tile grid: tile size, full tiles per dimension, tail of the last tile, tile count
static_assert(Solver::tile_size == 4 && Solver::full_tiles == 1 &&
Solver::last_tile_tail == 2 && Solver::num_tiles == 2);
// the knobs, read back from the solver type
static_assert(Solver::schedule == tdls::Schedule::LeftLooking);
static_assert(Solver::oot_threshold == 1e-10);
static_assert(Solver::singular_floor == std::numeric_limits<double>::min());
int piv[6];
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);