Compile-time dimension

TiledLUppSolverStatic on a 4 x 4 system in tiles of 2. The matrix, the right-hand sides and the expected solutions are declared above the shown code.

Factorize, then substitute

One factorization, two right-hand sides. The matrix holds the factors afterwards and the pivot array the row permutation.

\[\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_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}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];
double x1[4], x2[4];

// A becomes L and U in place, piv the row permutation; false means a singular matrix
if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// two right-hand sides on the same factorization
Solver::substitute<true, true, true>(A, 1, piv, 1, b1, x1, 1);
Solver::substitute<true, true, true>(A, 1, piv, 1, b2, x2, 1);

Factorize, then substitute in place

One buffer for the right-hand side and the 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 config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// y holds b on entry and x on exit: one buffer, permuted in place
Solver::substitute_inplace<true, true, true>(A, 1, piv, 1, y, 1);

Solve in one call

Factorization and substitution in one call, two buffers.

\[\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 config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];
double x[4];

// factorize and substitute in one call; A holds the factors afterwards
const bool ok = Solver::solve<true, true, true>(A, 1, piv, 1, b, x, 1);

Solve in place

Factorization and substitution in one call, one buffer. The forward substitution is folded into the factorization, so the result is bitwise the one of factorize then substitute in place.

\[\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 config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

// y holds b on entry and x on exit; the forward substitution runs inside the factorization
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);

One canonical column

The right-hand side is a canonical vector, generated on the fly: the solution is a column of the inverse.

\[\begin{split} A x = e_2, \quad 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 e_2 = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];
double x[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// x is column 2 of A^-1
Solver::substitute_canonical<true, true, true>(A, 1, piv, 1, 2, x, 1);

A block of canonical columns

Consecutive canonical columns solved together: the columns of the inverse a consistent tangent operator needs.

\[\begin{split} A X = \begin{pmatrix} e_0 & e_1 & e_2 \end{pmatrix}, \quad 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) \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// the first three columns of A^-1 in one sweep: column w of X starts at X + w * 4
double X[3 * 4];
Solver::substitute_canonical_multirhs<3, true, true, true>(A, 1, piv, 1, 0, X, 1, 0);

Factorize, then substitute a block

Three right-hand sides in one sweep, every tile loaded once for the block.

\[\begin{split} A X = B, \quad 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 & 7 & 8 \\ 14 & 7 & 3 \\ 24 & 8 & 6 \\ 33 & 10 & 5 \end{pmatrix}, \quad X = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// three right-hand sides in one sweep: column w of B and X starts at B + w * 4, X + w * 4
double X[3 * 4];
Solver::substitute_multirhs<3, true, true, true>(A, 1, piv, 1, B, X, 1, 0);

Factorize, then substitute a block in place

The same block, one buffer.

\[\begin{split} A X = B, \quad 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 & 7 & 8 \\ 14 & 7 & 3 \\ 24 & 8 & 6 \\ 33 & 10 & 5 \end{pmatrix}, \quad X = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// Y holds the three right-hand sides on entry and the three solutions on exit
Solver::substitute_inplace_multirhs<3, true, true, true>(A, 1, piv, 1, Y, 1, 0);

Solve a block

Factorization and block substitution in one call, two buffers.

\[\begin{split} A X = B, \quad 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 & 7 & 8 \\ 14 & 7 & 3 \\ 24 & 8 & 6 \\ 33 & 10 & 5 \end{pmatrix}, \quad X = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];
double X[3 * 4];

// factorize, then the three columns in one sweep
const bool ok = Solver::solve_multirhs<3, true, true, true>(A, 1, piv, 1, B, X, 1, 0);

Solve a block in place

Factorization and block substitution in one call, one buffer. The forward pass is not folded into the factorization here.

\[\begin{split} A X = B, \quad 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 & 7 & 8 \\ 14 & 7 & 3 \\ 24 & 8 & 6 \\ 33 & 10 & 5 \end{pmatrix}, \quad X = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 1 & 0 \\ 3 & 1 & 1 \\ 4 & 1 & 0 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

// Y holds the three right-hand sides on entry and the three solutions on exit
const bool ok = Solver::solve_inplace_multirhs<3, true, true, true>(A, 1, piv, 1, Y, 1, 0);

Pass cutting

A wide block cut into passes of two columns: a working-set knob, with the last pass taking the remainder.

\[\begin{split} A X = B, \quad B = \left(\begin{array}{cc|cc|c} 14 & 7 & 8 & 3 & 4 \\ 14 & 7 & 3 & 5 & 1 \\ 24 & 8 & 6 & 2 & 0 \\ 33 & 10 & 5 & 7 & 2 \end{array}\right), \quad X = \left(\begin{array}{cc|cc|c} 1 & 1 & 2 & 0 & 1 \\ 2 & 1 & 0 & 1 & 0 \\ 3 & 1 & 1 & 0 & 0 \\ 4 & 1 & 0 & 1 & 0 \end{array}\right) \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

if (!Solver::factorize<true, true>(A, 1, piv, 1)) return 1;

// five columns in passes of two, then the remaining one
double X[5 * 4];
Solver::substitute_multirhs<5, true, true, true, 2>(A, 1, piv, 1, B, X, 1, 0);

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} \textcolor{red}{10^{-12}} & 1 & 0 & 2 \\ \textcolor{red}{2 \cdot 10^{-12}} & 5 & 1 & 0 \\ \hline \textcolor{blue}{3} & 1 & 6 & 1 \\ 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}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

// one more argument: the columns whose best in-tile pivot fell below the threshold
int oot_count;
if (!Solver::factorize<true, true>(A, 1, piv, 1, oot_count)) return 1;

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

Singular verdict

Column 2, in red, is zero: no pivot reaches the floor and the call returns false.

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

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

int piv[4];

// column 2 is zero: false, and y is left partially updated
const bool ok = Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1);

Constant evaluation

The whole solve runs during constant evaluation. The data lives inside the constexpr function, so the whole program is shown.

\[\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 double solve_at_compile_time() {
    constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
    using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

    double A[4 * 4] = {4, 1, 0, 2, 1, 5, 1, 0, 0, 1, 6, 1, 2, 0, 1, 7};
    double y[4]     = {14, 14, 24, 33};
    int piv[4];

    if (!Solver::solve_inplace<true, true, true>(A, 1, piv, 1, y, 1)) return 0;
    return y[3];
}

// the whole solve runs during constant evaluation
constexpr double x3 = solve_at_compile_time();
static_assert(x3 > 4 - 1e-12 && x3 < 4 + 1e-12);

Newton iteration, then tangent columns

The MFront pattern on a small nonlinear system: one solve in place per iteration on a fresh jacobian, then the tangent columns from one factorization of the converged jacobian. The residual and jacobian functions are declared above the shown code.

\[\begin{split} F(x) = A x + x \circ x - c, \quad J(x) = A + 2\,\mathrm{diag}(x), \quad 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 c = \begin{pmatrix} 15 \\ 18 \\ 33 \\ 49 \end{pmatrix}, \quad x^\star = \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} \end{split}\]

constexpr auto config = tdls::TiledLUppConfig<double>{.tile_size = 2};
using Solver          = tdls::TiledLUppSolverStatic<double, 4, config>;

double x[4] = {0, 0, 0, 0};
double J[4 * 4], r[4];
int piv[4];

for (int iteration = 0; iteration < 20; ++iteration) {
    residual(A, c, x, r);
    jacobian(A, x, J);

    // one solve in place per iteration, on a fresh jacobian: r becomes the Newton step
    if (!Solver::solve_inplace<true, true, true>(J, 1, piv, 1, r, 1)) return 1;

    double step = 0;
    for (int i = 0; i < 4; ++i) {
        x[i] -= r[i];
        step = std::max(step, std::fabs(r[i]));
    }
    if (step < 1e-14) break;
}

// tangent columns: one factorization at the solution, two canonical columns in one block
jacobian(A, x, J);
if (!Solver::factorize<true, true>(J, 1, piv, 1)) return 1;

double dx_dc[2 * 4];
Solver::substitute_canonical_multirhs<2, true, true, true>(J, 1, piv, 1, 0, dx_dc, 1, 0);