Runtime dimension
TiledLUppSolverDynamic on a 5 x 5 system in tiles of 2: two full
tiles and a trailing tile of 1. The dimension, the matrix, the
right-hand sides and the expected solutions are declared above the
shown code, the arrays as std::vector<double>.
Factorize, then substitute
One factorization, two right-hand sides. The dimension is the first argument of every call, and every operand is a pointer plus a stride.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad
b_1 = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad
b_2 = \begin{pmatrix} 7 \\ 8 \\ 9 \\ 10 \\ 11 \end{pmatrix}, \quad
x_1 = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}, \quad
x_2 = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
std::vector<double> x1(n), x2(n);
// n is the first argument, every operand a pointer plus a stride; false means singular
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// two right-hand sides on the same factorization
Solver::substitute(n, A.data(), 1, piv.data(), 1, b1.data(), x1.data(), 1);
Solver::substitute(n, A.data(), 1, piv.data(), 1, b2.data(), x2.data(), 1);
Factorize, then substitute in place
One buffer for the right-hand side and the solution.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad b = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// y holds b on entry and x on exit
Solver::substitute_inplace(n, A.data(), 1, piv.data(), 1, y.data(), 1);
Solve in one call
Factorization and substitution in one call, two buffers.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad b = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
std::vector<double> x(n);
const bool ok = Solver::solve(n, A.data(), 1, piv.data(), 1, b.data(), x.data(), 1);
Solve in place
Factorization and substitution in one call, one buffer, the forward substitution folded into the factorization.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad b = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
// y holds b on entry and x on exit; the forward substitution runs inside the factorization
const bool ok = Solver::solve_inplace(n, A.data(), 1, piv.data(), 1, y.data(), 1);
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}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad e_2 = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
std::vector<double> x(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// x is column 2 of A^-1
Solver::substitute_canonical(n, A.data(), 1, piv.data(), 1, 2, x.data(), 1);
A block of canonical columns
Consecutive canonical columns solved together. The column count is a runtime argument, right after the dimension.
\[\begin{split}
A X = \begin{pmatrix} e_0 & e_1 & e_2 \end{pmatrix}, \quad A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right)
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// column w of X starts at X + w * n: element stride 1, column stride n
const int nrhs = 3;
std::vector<double> X(nrhs * n);
Solver::substitute_canonical_multirhs(n, nrhs, A.data(), 1, piv.data(), 1, 0, X.data(), 1, n);
Factorize, then substitute a block
Three right-hand sides in one sweep. A block carries two strides: the element stride within a column, and the stride between columns.
\[\begin{split}
A X = B, \quad A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad B = \begin{pmatrix} 12 & 7 & 6 \\ 16 & 8 & 2 \\ 27 & 9 & 7 \\ 40 & 10 & 2 \\ 50 & 11 & 10 \end{pmatrix}, \quad
X = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \\ 5 & 1 & 1 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// column w of B and X starts at B + w * n, X + w * n: element stride 1, column stride n
const int nrhs = 3;
std::vector<double> X(nrhs * n);
Solver::substitute_multirhs(n, nrhs, A.data(), 1, piv.data(), 1, B.data(), X.data(), 1, n);
Factorize, then substitute a block in place
The same block, one buffer.
\[\begin{split}
A X = B, \quad A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad B = \begin{pmatrix} 12 & 7 & 6 \\ 16 & 8 & 2 \\ 27 & 9 & 7 \\ 40 & 10 & 2 \\ 50 & 11 & 10 \end{pmatrix}, \quad
X = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \\ 5 & 1 & 1 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// Y holds the three right-hand sides on entry and the three solutions on exit
const int nrhs = 3;
Solver::substitute_inplace_multirhs(n, nrhs, A.data(), 1, piv.data(), 1, Y.data(), 1, n);
Solve a block
Factorization and block substitution in one call, two buffers.
\[\begin{split}
A X = B, \quad A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad B = \begin{pmatrix} 12 & 7 & 6 \\ 16 & 8 & 2 \\ 27 & 9 & 7 \\ 40 & 10 & 2 \\ 50 & 11 & 10 \end{pmatrix}, \quad
X = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \\ 5 & 1 & 1 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
const int nrhs = 3;
std::vector<int> piv(n);
std::vector<double> X(nrhs * n);
// factorize, then the three columns in one sweep
const bool ok =
Solver::solve_multirhs(n, nrhs, A.data(), 1, piv.data(), 1, B.data(), X.data(), 1, n);
Solve a block in place
Factorization and block substitution in one call, one buffer.
\[\begin{split}
A X = B, \quad A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad B = \begin{pmatrix} 12 & 7 & 6 \\ 16 & 8 & 2 \\ 27 & 9 & 7 \\ 40 & 10 & 2 \\ 50 & 11 & 10 \end{pmatrix}, \quad
X = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \\ 5 & 1 & 1 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
const int nrhs = 3;
std::vector<int> piv(n);
// Y holds the three right-hand sides on entry and the three solutions on exit
const bool ok =
Solver::solve_inplace_multirhs(n, nrhs, A.data(), 1, piv.data(), 1, Y.data(), 1, n);
Pass cutting
A wide block cut into passes of two columns. The column count stays a runtime argument, the pass width is a template argument.
\[\begin{split}
A X = B, \quad
B = \left(\begin{array}{cc|cc|c}
12 & 7 & 6 & 1 & 10 \\
16 & 8 & 2 & 6 & 2 \\
27 & 9 & 7 & 2 & 0 \\
40 & 10 & 2 & 8 & 0 \\
50 & 11 & 10 & 1 & 2
\end{array}\right), \quad
X = \left(\begin{array}{cc|cc|c}
1 & 1 & 1 & 0 & 2 \\
2 & 1 & 0 & 1 & 0 \\
3 & 1 & 1 & 0 & 0 \\
4 & 1 & 0 & 1 & 0 \\
5 & 1 & 1 & 0 & 0
\end{array}\right)
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1)) return 1;
// passes of two: nrhs at run time, the pass width as template argument
const int nrhs = 5;
std::vector<double> X(nrhs * n);
Solver::substitute_multirhs<2>(n, nrhs, A.data(), 1, piv.data(), 1, B.data(), X.data(), 1, n);
Out-of-tile counter
The two candidates of column 0 inside the first tile, in red, are below the threshold. The search below the tile takes row 2, in blue, and the counter reports the column.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
\textcolor{red}{10^{-12}} & 1 & 0 & 0 & 1 \\
\textcolor{red}{2 \cdot 10^{-12}} & 6 & 1 & 0 & 0 \\ \hline
\textcolor{blue}{3} & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad
b = \begin{pmatrix} 7 \\ 15 \\ 27 \\ 40 \\ 49 \end{pmatrix}, \quad
x = \begin{pmatrix} 0 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
std::vector<int> piv(n);
// one more argument: the columns whose best in-tile pivot fell below the threshold
int oot_count;
if (!Solver::factorize(n, A.data(), 1, piv.data(), 1, oot_count)) return 1;
const bool searched_below_the_tile = oot_count == 1 && piv[0] == 2;
Tile helpers
The grid of a runtime dimension, read from the solver: the number of tiles per dimension and the size of the tile at a position.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad b = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver = tdls::TiledLUppSolverDynamic<double, config>;
// the grid of a runtime dimension: three tiles per dimension for n = 5, the last of size 1
const int tiles = Solver::num_tiles(n);
const int last = Solver::tile_size_at((tiles - 1) * Solver::tile_size, n);
std::vector<int> piv(n);
const bool ok = Solver::solve_inplace(n, A.data(), 1, piv.data(), 1, y.data(), 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::TiledLUppSolverDynamic<double, wide>;
std::vector<int> piv(n);
const bool ok = Solver::solve_inplace(n, A.data(), 1, piv.data(), 1, y.data(), 1);
Same results as the compile-time solver
The two solvers at equal shape and configuration: the factors, the pivots and the solution are bitwise identical.
\[\begin{split}
A = \left(\begin{array}{cc|cc|c}
5 & 1 & 0 & 0 & 1 \\
1 & 6 & 1 & 0 & 0 \\ \hline
0 & 1 & 7 & 1 & 0 \\
0 & 0 & 1 & 8 & 1 \\ \hline
1 & 0 & 0 & 1 & 9
\end{array}\right), \quad b = \begin{pmatrix} 12 \\ 16 \\ 27 \\ 40 \\ 50 \end{pmatrix}, \quad x = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{pmatrix}
\end{split}\]
constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Dynamic = tdls::TiledLUppSolverDynamic<double, config>;
using Static = tdls::TiledLUppSolverStatic<double, 5, config>;
std::vector<int> piv_dynamic(n);
int piv_static[5];
// same shape, same configuration: the same arithmetic, bit for bit
const bool ok_dynamic =
Dynamic::solve_inplace(n, A_dynamic.data(), 1, piv_dynamic.data(), 1, y_dynamic.data(), 1);
const bool ok_static =
Static::solve_inplace<true, true, true>(A_static, 1, piv_static, 1, y_static, 1);