decomposition.hpp File Reference

Matrix factorisations and the algorithms built on them. More...

#include <algorithm>
#include <array>
#include <atomic>
#include <bit>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <functional>
#include <limits>
#include <memory>
#include <numbers>
#include <optional>
#include <stdexcept>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
#include "matrix.hpp"
#include "simd.hpp"
#include "parallel.hpp"
#include "vector.hpp"

Include dependency graph for decomposition.hpp:

Go to the source code of this file.

Classes
class	numerical::singular_matrix_error
	Thrown when a matrix has no unique solution (it is singular). More...
class	numerical::non_finite_error
	Thrown when a decomposition's input holds an infinite or NaN entry: no factorisation of it means anything, and the iterative ones (eigh, SVD) would stop at once and return plausible values. More...
class	numerical::lu_decomposition< T >
	An LU factorisation with partial pivoting of a square matrix. More...
class	numerical::not_positive_definite_error
	Thrown when a Cholesky factorisation finds the matrix not positive-definite. More...
class	numerical::cholesky_decomposition< T >
	A Cholesky factorisation A = L Lᴴ of a Hermitian positive-definite matrix. More...
struct	numerical::detail::thin_qr< T >
	A thin QR factorisation, row-major: q is rows × cols with orthonormal columns, r is cols × cols upper triangular with a real, non-negative diagonal. More...
struct	numerical::detail::reflector_block< T >
	The reflectors of one panel of the blocked QR, I − V T Vᴴ as LAPACK's dlarft forms it. More...
struct	numerical::detail::qr_reflections< T >
	The reflections of a Householder QR, Q left as its blocks of reflectors. More...
class	numerical::rank_cutoff
	The magnitude at or under which a singular value counts as zero, in a rank (matrix_rank), a pseudo-inverse (pinv) or a least-squares solution (lstsq): by default NumPy's, the largest singular value times the machine epsilon times the larger dimension; relative to the largest, as NumPy's rcond and rtol; or absolute, as matrix_rank's tol. More...
class	numerical::qr_decomposition< T >
	A QR factorisation A = Q R of an m × n matrix, as NumPy's qr, reduced by default. More...
class	numerical::convergence_error
	Thrown when an iterative algorithm has not converged within its bound of sweeps. More...
struct	numerical::detail::tridiagonal_reduction< T >
	The unitary M = Q Dᴴ of a reduction to real tridiagonal form, A = M T Mᴴ, kept as its factors. More...
struct	numerical::detail::row_operation
	One operation of tridiagonal_ql or bidiagonal_qr on the rows of a vector matrix, recorded to be applied later. More...
struct	numerical::detail::rotation_run
	A run of recorded rotations in which each shares a row with the next, as a sweep of the QL or QR iteration gives them: rows (i, i+1) for i going down (down), or going up. More...
struct	numerical::detail::secular_solution
	A root of a secular equation: the pole it was sought from, and its offset from that pole. More...
struct	numerical::detail::workspace_buffer
	One buffer of merge_workspace: its storage and how many elements it holds. More...
struct	numerical::detail::merge_workspace
	The buffers of merge_rank_one, kept from one join to the next: a fresh buffer per join is mapped anew, page by page. Each is written before it is read (updated zeroed row by row by the gathers). More...
struct	numerical::detail::divide_node
	One block of divide_and_conquer's tree: its rows [0, n) of the diagonal d its parent reads, the diagonal and subdiagonal it is solved on (its own copies scaled by a power of two, or d and the parent's), and its diagonal block of the eigenvector matrix, q, rows ldq apart. More...
class	numerical::symmetric_eigen
	The eigenvalues and eigenvectors of a real symmetric matrix. More...
class	numerical::hermitian_eigen
	The eigenvalues and eigenvectors of a complex Hermitian matrix. More...
struct	numerical::detail::arrow_workspace
	The buffers merge_arrow works in, kept by each thread from one join to the next (thread_arrow_workspace): mapping their pages anew took 11–16 % of an svd of order 64 to 512 on one x86-64 core, 13–20 % on six. More...
struct	numerical::detail::bidiagonal_node
	One block of bidiagonal_divide's tree: the bidiagonal it decomposes (its own copies scaled by a power of two, or its parent's), the two halves it splits into, and its decomposition once solved. More...
class	numerical::svd_decomposition< T >
	The singular value decomposition A = U Σ Vᴴ of an m × n matrix. More...
struct	numerical::least_squares_solution< T >
	A least-squares solution and what numpy.linalg.lstsq reports beside it (least_squares). More...

Namespaces
namespace	numerical
	The NumericAL public API.
namespace	numerical::detail

Typedefs
template<typename T >
using	numerical::detail::magnitude_t = decltype(std::abs(std::declval< const T & >()))
	The type std::abs gives for a T: the real type of a complex one, T itself otherwise.
template<typename T >
using	numerical::detail::value_range_t = std::conditional_t< std::is_same_v< T, float >||std::is_same_v< T, std::complex< float > >, float, double >
	The real type whose range a T's values take: float for a single-precision scalar, real or complex; double otherwise (a long double's magnitudes are carried in double).
using	numerical::detail::symmetric_row_kernel = double(*)(const double *, const double *, double, double *, std::size_t)
	A kernel of symmetric_row's form: symmetric_row itself, or an x86-64 one.
using	numerical::detail::symmetric_rows_kernel = std::array< double, symmetric_row_group >(*)(const double *, std::size_t, const double *, double *, std::size_t)
	A kernel of symmetric_rows' form: symmetric_rows itself, or an x86-64 one.

Enumerations
enum class	numerical::qr_mode { numerical::reduced , numerical::complete }
	The factors a numerical::qr_decomposition keeps, as numpy.linalg.qr's modes name them. More...
enum class	numerical::eigen_method { numerical::tridiagonal_ql , numerical::jacobi , numerical::divide_and_conquer }
	How numerical::symmetric_eigen diagonalises. More...
enum class	numerical::svd_method { numerical::bidiagonal_qr , numerical::jacobi , numerical::divide_and_conquer }
	How numerical::svd_decomposition computes the decomposition of R, after a Householder QR. More...

Functions
template<typename T >
T	numerical::detail::pivot_divisor (const T &pivot)
	The divisor over_pivot takes for the pivot pivot. A complex pivot's reciprocal, once, as LAPACK's zgetf2 takes it (a complex division calls the C library's robust one, a quarter of a complex LU of order 48 when it ran for every element); a real pivot itself, each element divided.
template<typename R >
std::complex< R >	numerical::detail::pivot_divisor (const std::complex< R > &pivot)
	pivot_divisor of a complex pivot: its reciprocal.
template<typename T >
T	numerical::detail::over_pivot (const T &x, const T &divisor)
	x over the pivot whose pivot_divisor is divisor: the same quotient wherever the LU takes it.
template<typename R >
std::complex< R >	numerical::detail::over_pivot (const std::complex< R > &x, const std::complex< R > &divisor)
	over_pivot of a complex x: times the pivot's reciprocal, the product in line.
template<typename T >
void	numerical::detail::subtract_term (T &sum, const T &a, const T &b)
	sum −= a · b in one statement, which the compiler may fuse (as it fuses the matrix solve's).
template<typename R >
void	numerical::detail::subtract_term (std::complex< R > &sum, const std::complex< R > &a, const std::complex< R > &b)
	subtract_term of complex numbers: the product in line (operator_times, the operator's rounding).
constexpr std::size_t	numerical::detail::lu_block_for (std::size_t n) noexcept
	The panel width of the LU of order n: a power of two, so a row's panel is a mask away, not a division.
constexpr std::size_t	numerical::detail::cholesky_block_for (std::size_t n) noexcept
	The panel width of the Cholesky factorisation of order n: narrower up to order 256, where the panels, row by row on one thread, are most of the time; a power of two, which the product kernels' tiles divide (widths of 24, 40 or 56 left the AMX's 32 × 16 tiles part empty: order 320, 20 % longer).
template<typename T >
T	numerical::detail::conj_of (const T &x)
	The conjugate of a real (or general) scalar — the value itself.
template<typename U >
std::complex< U >	numerical::detail::conj_of (const std::complex< U > &x)
	The conjugate of a complex scalar.
template<typename T >
T	numerical::detail::from_double (double x)
	x, a magnitude carried in double, as the scalar T: converted to T's real type first, so a std::complex<float> is built from a float (a double argument narrows inside its constructor).
template<typename T >
double	numerical::detail::real_part (const T &x)
	The real part of a real (or general) scalar as a double.
template<typename U >
double	numerical::detail::real_part (const std::complex< U > &x)
	The real part of a complex scalar.
template<typename T >
double	numerical::detail::rescue_factor (double largest)
	The exact power of two that brings a magnitude below about the square root of the smallest normal number of T's range (2⁻⁵⁰⁰ in double, 2⁻⁶⁰ in float) up among normal numbers, as LAPACK's dlarfg rescales; 1 for any larger magnitude. It is a T itself: a zero takes the factor of the smallest normal number, which a double's range holds and a float's does not.
template<typename T >
T	numerical::detail::phase_of (const T &x)
	The unit-modulus factor of x: x / |x| (its sign, for a real scalar), and 1 at 0. A subnormal x is scaled up first, exactly: its modulus would round to the coarse subnormal grid, and x / |x| miss unit modulus by as much.
template<typename T >
double	numerical::detail::rescale_tiny (buffer< T > &x)
	Scales x by rescue_factor of its largest magnitude, exactly, so that a Householder vector built from it is computed among normal numbers: a subnormal norm rounds to the coarse subnormal grid, and the reflector is then not orthogonal.
template<typename T >
double	numerical::detail::norm_at_scale (const T *x, std::size_t count)
	The Euclidean norm of the count elements at x, summed at the scale of their largest magnitude: a plain sum of squares overflows near 1e155 and underflows to 0 near 1e-155.
template<typename T >
double	numerical::detail::plain_or_at_scale (double plain, const T *x, std::size_t count)
	√plain, the plain sum of the squares of the count elements at x, where it neither overflowed nor underflowed (between plain_squares_low and plain_squares_high); else norm_at_scale.
template<typename T >
double	numerical::detail::scaled_norm (const T *x, std::size_t count)
	The Euclidean norm of the count elements at x, by norm_at_scale for this type.
double	numerical::detail::scaled_norm (const double *x, std::size_t count)
	The same for double: the plain sum of squares in one pass of dot_product where it holds.
double	numerical::detail::scaled_norm (const std::complex< double > *x, std::size_t count)
	The same for std::complex<double>: dot_product over the parts, which an array of them holds in order ([complex.numbers]).
template<typename T >
double	numerical::detail::scaled_norm (const buffer< T > &x)
	scaled_norm of a whole vector.
double	numerical::detail::rotation_length (double a, double b) noexcept
	√(a² + b²) for a plane rotation of the QR and QL iterations, as LAPACK's dlartg takes it: the plain sum of squares where the larger magnitude lies within 2^±500, so its square neither overflows nor underflows and the smaller one's underflow is below its rounding; std::hypot, which scales every time, outside. A graded matrix reaches both ends though scaled to its largest entry: a QL shift's ratio (d₁ − d₀) / 2e came to 5e199 on [[1e-200, 1e-200], [1e-200, 1]].
template<typename T >
bool	numerical::detail::is_finite_entry (const T &x)
	Whether x is finite.
template<typename U >
bool	numerical::detail::is_finite_entry (const std::complex< U > &x)
	Whether x is finite: both of its parts.
template<typename T >
void	numerical::detail::require_finite (const matrix< T > &a, const char *what)
	The one non-finite policy, checked where every decomposition starts.
bool	numerical::detail::all_finite (const double *x, std::size_t count) noexcept
	Whether the count values at x are all finite: a NaN or an infinity has every exponent bit set, the one exponent that one more carries into the sign bit. One pass with no early exit and no comparison, an integer AND, add and OR (a check per element and an exit at the first was a tenth of an LU of order 128; the comparison per element, one M1 Pro core: 6.4 → 8.5 values a nanosecond). Under GCC and Clang eight values a step in four vectors of two, whatever the compiler would vectorise: GCC at -O2 (the arm64 wheels') vectorised the plain loop not at all, at -O3 by one vector (1089 values, one arm64 core: GCC -O2 541 → 132 ns, -O3 299 → 132, Clang level at 132; one x86-64 core, GCC 13 -O3 271 → 163 ns, Clang 18 level); a Cholesky factorisation of order 33 lost a fifth of its time to the plain loop under GCC -O2.
template<typename T >
matrix< T >	numerical::detail::parallel_copy_of (const matrix< T > &a)
	A copy of a whose pages the threads fault in (parallel_copy), a decomposition's own storage: the matrix's copy faulted every one in on the calling thread (order 1000 on one arm64 core 0.64 ms, an eighth of the Cholesky factorisation; on eight cores 0.28 ms).
bool	numerical::detail::all_finite_shared (const double *x, std::size_t count)
	all_finite on the threads by chunks of finite_chunk, past two of them: a pass bound by the memory, an element read weighed as half one copied (six x86-64 AVX2 cores, a Cholesky factorisation of order 500 with its other passes on the threads too: 1.12 → 0.91 ms).
void	numerical::detail::require_finite (const matrix< double > &a, const char *what)
	require_finite for double, its elements contiguous: all_finite.
void	numerical::detail::require_finite (const matrix< std::complex< double > > &a, const char *what)
	The same for std::complex<double>, whose array holds its parts in order ([complex.numbers]).
template<typename T >
bool	numerical::detail::lower_is_finite (const matrix< T > &a)
	Whether the lower triangle of the square a, and its diagonal's real part, are finite: entry by entry, kept out of line: the rare road of require_finite_lower where the whole matrix is not finite.
template<typename T >
void	numerical::detail::require_finite_lower (const matrix< T > &a, const char *what)
	require_finite over what a Hermitian eigensolver reads: the lower triangle, and the diagonal's real part (as numpy.linalg.eigh reads them, whatever the rest holds).
void	numerical::detail::require_finite_lower (const matrix< double > &a, const char *what)
	require_finite_lower for double: all_finite over the whole matrix, the one vectorised pass require_finite takes, and the triangle entry by entry only where it finds a NaN or an infinity. Out of line, so that a decomposition's constructor stays small enough to be inlined: in it, the check kept GCC from inlining the Cholesky's (order 8, 0.204 → 0.220 µs), as an entry-by-entry check (6 % at order 64 on arm64) or a row-by-row one (21 % at order 8 under GCC) cost in time.
void	numerical::detail::require_finite_lower (const matrix< std::complex< double > > &a, const char *what)
	The same for std::complex<double>, whose array holds its parts in order ([complex.numbers]).
template<typename T >
void	numerical::detail::subtract_panel_columns (T *column, const T *columns, std::size_t height, const T *l, std::size_t terms)
	column −= Σₜ x_t · conj(l[t]) over the terms columns columns (height elements each, one after another), term by term in ascending order: an elimination step of a Cholesky panel's rows below (an overload per scalar, as conj_of). Each term subtract_product, the column read and written again for every one.
void	numerical::detail::subtract_panel_columns (double *column, const double *columns, std::size_t height, const double *l, std::size_t terms)
	subtract_panel_columns for double: subtract_products, which holds the column in registers through the terms, the same multiply-adds in the same order, so the same values (a complex one keeps the term-by-term calls, whose x · f the other would round as f · x).
template<typename T >
vector< T >	numerical::solve (matrix< T > a, const vector< T > &b)
	Solves the linear system A x = b for a square A.
template<typename T >
matrix< T >	numerical::solve (matrix< T > a, const matrix< T > &b)
	Solves A X = B for a square A and every column of B at once, as numpy.linalg.solve with a matrix right-hand side: one factorisation of A, then one substitution per column.
template<typename T >
T	numerical::determinant (matrix< T > a)
	The determinant of a square matrix.
template<typename T >
matrix< T >	numerical::inverse (matrix< T > a)
	The inverse of a square matrix.
template<typename T >
cholesky_decomposition< T >	numerical::cholesky (matrix< T > a)
	Cholesky-factorises a Hermitian positive-definite matrix.
template<typename T >
void	numerical::detail::reflect_rows (buffer< T > &a, std::size_t cols, const buffer< T > &v, std::size_t first, std::size_t from)
	Applies I − 2 v vᴴ from the left to the rows first…first+|v|-1 of a row-major array, on its columns from…cols-1.
template<typename T >
reflector_block< T >	numerical::detail::gather_reflectors (const std::vector< buffer< T > > &reflectors, std::size_t k0, std::size_t k1, std::size_t rows, const std::vector< T > *taus=nullptr)
	Gathers the reflectors reflectors[k0…k1-1] (each I − τ_k v vᴴ: τ_k from taus, else 2 for a unit v; an empty one the identity) into I − V T Vᴴ = H_{k0} H_{k0+1} … H_{k1-1}, T as LAPACK's dlarft forms it.
template<typename T >
void	numerical::detail::apply_reflector_block (const reflector_block< T > &block, bool adjoint, buffer< T > &a, std::size_t cols, std::size_t from, std::size_t to)
	Applies I − V T Vᴴ (or its adjoint I − V Tᴴ Vᴴ) from the left to the block's rows of a row-major array, on its columns from…to-1: W = Vᴴ A, then W ← T W or Tᴴ W, then A −= V W, the two outer products through the matrix product kernel.
template<typename T >
void	numerical::detail::apply_transposed_reflector_block (const reflector_block< T > &block, buffer< T > &at, std::size_t cols, std::size_t from, std::size_t to, bool adjoint=false)
	apply_reflector_block to A held transposed: at is Aᵀ (row-major, cols columns), and A ← (I − V T Vᴴ) A on A's rows first… and columns from…to-1 is Aᵀ ← Aᵀ (I − conj(V) Tᵀ Vᵀ) on Aᵀ's columns first… and rows from…to-1 (with the adjoint (I − V Tᴴ Vᴴ) A, conj(T) for Tᵀ). Its three products have to − from rows, which the matrix product kernel runs well and on every core; apply_reflector_block's first has the block's few (a reflector block's 32 rows by 1000 deep: 54 GFLOP/s on twenty x86-64 threads).
template<typename T > requires std::is_floating_point_v<T>
T	numerical::detail::conjugate_dot (const T *x, const T *y, std::size_t count)
	Σ conj(x[i]) · y[i] for a real T, where the conjugate is the value: dot_product.
template<typename T >
T	numerical::detail::conjugate_dot (const T *x, const T *y, std::size_t count)
	Σ conj(x[i]) · y[i] of complex numbers (conjugate_dot_product).
template<typename T >
void	numerical::detail::reflect_columns (T *panel, std::size_t height, const buffer< T > &v, std::size_t first, std::size_t from, std::size_t to, const T &s=T {2})
	Applies I − s v vᴴ from the left (s 2 for a unit v) to the columns from…to-1 of the column-major panel (height rows per column), on its rows first…first+|v|-1: each column a dot product with v, then v scaled out of it, both along the column at unit stride.
template<typename T >
void	numerical::detail::factor_qr_panel (T *panel, std::size_t ld, std::size_t top, std::vector< buffer< T > > &reflectors, std::vector< T > &taus, std::size_t k0, std::size_t k1)
	Householder reflections of the columns k0…k1-1 of a column-major matrix, rows from k0, one by one, in place, into reflectors[k0…k1-1] and taus: panel holds column k0, ld elements per column from row top (0 in the whole matrix, k0 in a copy of the panel). A reflector meets each column at unit stride. Recursive halves, the left's reflectors applied to the right at once, took longer at these widths (32 columns).
template<typename T >
std::vector< reflector_block< T > >	numerical::detail::reflect_in_columns (buffer< T > &columns, std::size_t rows, std::size_t cols)
	The blocks of reflectors of a Householder QR of the column-major columns (rows per column, cols columns), factorised in place as LAPACK's dgeqrf leaves it: R on and above the diagonal, each reflector's w₁… below it. The columns right of a panel take its block as (I − V T Vᴴ)ᴴ, held transposed (apply_transposed_reflector_block): the products then have as many rows as there are columns to update, where the row-major form's first product has the block's 32 (one x86-64 core: 39 ms against 46 for the updates of order 1000).
template<typename T >
void	numerical::detail::transpose_into (const T *from, std::size_t rows, std::size_t cols, T *to)
	The row-major rows × cols from into to as its transpose, cols × rows row-major, by tiles, bands of them on the threads (a pass bound by the memory, its bounds captured by value).
template<typename T >
void	numerical::detail::factor_row_major_panel (buffer< T > &a, std::size_t rows, std::size_t cols, std::vector< buffer< T > > &reflectors, std::vector< T > &taus, std::size_t k0, std::size_t k1)
	factor_qr_panel on the columns k0…k1-1 of the row-major a, rows from k0, in a column-major copy of the panel, copied back. A function of its own: written into reflect_in_rows' loop, the same instructions took 5 % more cycles at order 128 (clang, arm64).
template<typename T >
std::vector< reflector_block< T > >	numerical::detail::reflect_in_rows (buffer< T > &a, std::size_t rows, std::size_t cols, buffer< T > *r=nullptr)
	The blocks of reflectors of a Householder QR of the row-major a (rows × cols), factorised in place as reflect_in_columns leaves the column-major one (R on and above the diagonal, each reflector's w₁… below it): each panel in a column-major copy, the columns right of it updated by apply_reflector_block, whose projection reads A as the product's right operand. The AMX tile reads that in place, and NEON's packs it faster than the column-major form's left one (arm64, order 512: 6.3 ms against 7.2 under clang, 16.4 against 17.9 under GCC). Where r is given, it is R (min(rows, cols) × cols row-major, zeros below the diagonal).
template<typename T >
qr_reflections< T >	numerical::detail::householder_reflect (buffer< T > a, std::size_t rows, std::size_t cols)
	The reflections of householder_qr, without forming Q: A = Q R with Q the product of the blocks' reflectors applied to the first k = min(rows, cols) columns of the identity, R's diagonal not yet made real (a product with Qᴴ, or with Q of a block padded with zeros, applies the blocks instead). a is row-major. On the x86-64 units it is transposed in by tiles and factorised column-major (reflect_in_columns, whose updates' products run faster there: one AVX2 core, order 1000, 39 ms against 46), elsewhere in place (reflect_in_rows).
template<typename T >
void	numerical::detail::form_q_into (buffer< T > &out, const std::vector< reflector_block< T > > &blocks, std::size_t rows, std::size_t first, std::size_t count, bool transposed)
	The columns [first, first + count) of Q = H₀ H₁ …, the blocks applied (the last first) to those of the identity, into out: rows × count row-major, or where transposed Qᵀ row-major (Q column-major). Formed as householder_reflect factorises, so that a Q is the same to the bit in either layout: as the rows of Qᵀ on the x86-64 units (apply_transposed_reflector_block), row-major elsewhere (apply_reflector_block), then transposed where the other layout is asked. The columns of the identity left of a block's rows are left as they are.
template<typename T >
thin_qr< T >	numerical::detail::householder_qr (buffer< T > a, std::size_t rows, std::size_t cols, bool complete=false)
	Factorises the rows × cols row-major a by Householder reflections, A = Q R with k = min(rows, cols): Q is rows × k with orthonormal columns and R k × cols upper trapezoidal. Each column below the diagonal is reflected onto it, so Q is orthonormal to rounding whatever the conditioning of A, and a column dependent on the earlier ones leaves a zero on R's diagonal.
double	numerical::detail::rank_tolerance (double largest, std::size_t rows, std::size_t cols)
	The magnitude under which a diagonal entry of R or a singular value counts as zero, as NumPy draws it: the largest times the machine epsilon times the larger dimension.
template<typename T >
qr_decomposition< T >	numerical::qr (const matrix< T > &a)
	QR-factorises an m × n matrix.
template<typename T > requires std::is_floating_point_v<T>
void	numerical::detail::add_conjugate_product (T *y, const T *x, const T &f, std::size_t count)
	y[i] += conj(x[i]) · f for a real T, where the conjugate is the value: add_product.
template<typename T >
void	numerical::detail::add_conjugate_product (T *y, const T *x, const T &f, std::size_t count)
	y[i] += conj(x[i]) · f, each product's statements spelled out as the complex operator's (four part products rounded apart, ac − bd and ad + bc, then their sums with y): the same bits for finite parts whose products stay finite, where the operator's rescue of a product whose parts are both NaN never acts, and a loop the compiler vectorises (the reduction's rows, scaled by their largest magnitude, are such).
template<typename T >
T	numerical::detail::hermitian_row (symmetric_row_kernel, const T *row, const T *v, T *out, std::size_t i)
	Row i's terms of a Hermitian product from the lower triangle: out[j] += conj(row[j]) v[i] for j < i, and the return Σⱼ≤ᵢ row[j] v[j], for a complex T.
std::complex< double >	numerical::detail::hermitian_row (symmetric_row_kernel kernel, const std::complex< double > *row, const std::complex< double > *v, std::complex< double > *out, std::size_t i)
	The same for a std::complex<double> row: in one pass on AVX-512 or AVX2 with FMA where the matrix product's kernel runs there (avx512_hermitian_row), else as the template.
double	numerical::detail::hermitian_row (symmetric_row_kernel kernel, const double *row, const double *v, double *out, std::size_t i)
	The same for a real row, in one pass through kernel.
template<typename T >
std::size_t	numerical::detail::symmetric_row_speedup ()
	The weight of a symmetric row's multiply-adds in hermitian_product: on the x86-64 units, where its rows are paced by one core's loads, vector_speedup, as more threads pay than their multiply-adds alone would pay for (an i9-9900K's reduction of order 512 took 17 % less, an i9-10900X's 7 %); 1 elsewhere (arm64's, weighed so, took as long or longer).
template<typename T >
void	numerical::detail::hermitian_product (const T *block, std::size_t n, std::size_t m, const T *v, T *w, buffer< T > &scratch)
	w = B v for the Hermitian m × m block (row stride n) from its lower triangle alone, as LAPACK's dsymv: each row's lower part is read once and serves twice, w[i] += Σⱼ≤ᵢ B[i][j] v[j] and w[j] += conj(B[i][j]) v[i] for j < i, so a reduction step reads half the block it would read whole. Rows go by tasks of fixed size, each into its own partial w, summed in a fixed order: the result does not depend on the number of threads. A real row is read once for both (symmetric_row), on the widest x86-64 units the matrix product's kernel may use (gemm_path_current).
template<typename T > requires std::is_floating_point_v<T>
void	numerical::detail::subtract_corrections (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_products of a real T, where the reduction's corrections are those of its type.
template<split_complex T>
void	numerical::detail::subtract_corrections (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_products of complex numbers, each term y[i] −= f · x[i] spelled out as the complex operator's statements (four part products rounded apart, ac − bd and ad + bc, then their differences with y): the same bits for finite parts whose products stay finite, where the operator's rescue of a product whose parts are both NaN never acts, and loops the compiler vectorises (the reduction's rows, scaled by their largest magnitude, are such). Complex double on the x86-64 units keeps its fused kernel.
template<typename T >
void	numerical::detail::corrected_hermitian_product (const T *block, std::size_t n, std::size_t m, const T *v, T *w, buffer< T > &scratch, const T *corrections, std::size_t ldc, std::size_t terms, T *factors, T *next=nullptr, const T *next_factors=nullptr)
	hermitian_product corrected by a panel's earlier steps: w = B v − Σₜ fₜ cₜ, the coefficients f[2p] = (c[2p+1] · v) and f[2p+1] = (c[2p] · v) (conjugated) computed beside the product, among its tasks, and the partials summed and the corrections taken by bands of rows: each element the operations, in the order, hermitian_product, conjugate_dot and subtract_products give it apart.
double	numerical::detail::real_subdiagonal (double t, double, double &)
	A real subdiagonal entry of the tridiagonal: kept with its sign (the rotations take any sign).
double	numerical::detail::real_subdiagonal (std::complex< double > t, std::complex< double > delta, std::complex< double > &next)
	A complex subdiagonal entry made real by the next row's unit phase: next = delta · conj(phase(t)).
template<typename T >
tridiagonal_reduction< T >	numerical::detail::tridiagonalize (buffer< T > &a, std::size_t n, buffer< double > &d, buffer< double > &e)
	Reduces the Hermitian (real: symmetric) n × n row-major a, of which the lower triangle is read, to a real symmetric tridiagonal by Householder reflections, A = M T Mᴴ with M unitary.
template<typename T >
buffer< T >	numerical::detail::tridiagonal_basis (const tridiagonal_reduction< T > &reduction, std::size_t n)
	Mᵀ of a reduction, row-major (its rows are the columns of M): Qᵀ by blocks of reflectors from the last back, as LAPACK's dorgtr forms Q, each block on the rows it reaches of the identity, then row j times conj(δⱼ).
template<typename T >
void	numerical::detail::to_matrix_basis (const tridiagonal_reduction< T > &reduction, buffer< T > &at, std::size_t n)
	The rows of at, vectors in the tridiagonal's basis, taken to the matrix's: at ← at Mᵀ = at Dᴴ Qᵀ, column j times conj(δⱼ), then the blocks of reflectors from the last back on every row, as LAPACK's dormtr: two thirds of the operations of forming M (tridiagonal_basis) and multiplying by it.
std::vector< rotation_run >	numerical::detail::rotation_runs (const std::vector< row_operation > &operations)
	The operations split into runs of chained rotations, every other operation a run of one.
template<typename T >
void	numerical::detail::apply_row_operations (buffer< T > &m, std::size_t n, const std::vector< row_operation > &operations)
	Applies the recorded operations, in order, to the row-major m (n per row). Each column takes them independently of the others, so chunks of columns run in parallel, every element receiving the same operations in the same order as applying each in place would give it.
template<typename T >
void	numerical::detail::tridiagonal_ql (buffer< double > &d, buffer< double > &e, buffer< T > &zt, std::size_t n, std::size_t max_steps)
	The eigenvalues of the symmetric tridiagonal (d, e) by implicit QL with Wilkinson-type shifts (EISPACK's tql2), the rotations accumulated into the rows of zt.
double	numerical::detail::quadratic_step (double a, double b, double c)
	The smaller root η of C η² − A η + B = 0, in the form that does not cancel: a step of solve_secular's fit (shared by its instantiations, whose data reach one form or the other).
template<typename Offset >
secular_solution	numerical::detail::solve_secular (const buffer< double > &z, double rho, std::size_t i, std::size_t max_iterations, double *delta, const Offset &offset)
	Root i of 1/ρ + Σⱼ zⱼ²/(Pⱼ − x) = 0, poles P strictly ascending, no zⱼ zero, ρ > 0, as LAPACK's dlaed4 and dlasd4 find it. The poles are given only as their offsets from one another, offset(o, j) = Pⱼ − P_o, which the caller computes accurately.
double	numerical::detail::secular_root (const buffer< double > &d, const buffer< double > &z, double rho, std::size_t i, std::size_t max_iterations, double *delta)
	Root i of 1/ρ + Σⱼ zⱼ²/(dⱼ − λ) = 0 (d strictly ascending, ‖z‖ = 1), the i-th eigenvalue of diag(d) + ρ z zᵀ, by solve_secular.
double *	numerical::detail::workspace_of (workspace_buffer &buffer, std::size_t size)
	buffer at least size elements long, its storage kept for the next join. A larger one is mapped fresh, and its pages are touched by every thread, chunk by chunk: on Linux a page costs about a microsecond the first time it is written, and the join that first wrote it waited on them one thread at a time (an eigh of order 1024 mapped some 30 MB of these buffers on one x86-64 core).
merge_workspace &	numerical::detail::thread_workspace ()
	The calling thread's merge_workspace, kept from one decomposition to the next: mapping its pages anew took 8–11 % of the eigenvectors' time (six AVX2 cores, an eigh of order 256 or 512). A join uses it whole on one thread, and nests no other join; release_large_workspace lets a large buffer go after.
void	numerical::detail::release_large_workspace (merge_workspace &work)
	Releases the buffers of work longer than kept_workspace.
void	numerical::detail::merge_rank_one (buffer< double > &d, buffer< double > z, double rho, double *q, std::size_t ldq, std::size_t n, std::size_t split, merge_workspace &work)
	The eigenvalues of diag(d) + ρ z zᵀ, ascending, into d, and its eigenvectors times q (the row-major n × n block of rows ldq apart whose columns span the halves) into q, as LAPACK's dlaed1.
void	numerical::detail::plan_divide (std::vector< divide_node > &nodes, double *d, const double *e, std::size_t n, double *q, std::size_t ldq, std::size_t depth)
	Appends the block (d, e, n) and, before it is solved, everything below it to nodes, parents before children: each scaled to work near unit size by a power of two (exactly), as the deflation tolerance is relative to the weights, which have unit norm, and so to the tridiagonal only at that size; each half has its own. A block past divide_leaf is split into two halves plus a rank-one coupling, as LAPACK's dlaed0: T = diag(T₁ − ρ e_a e_aᵀ, T₂ − ρ e_b e_bᵀ) + ρ u uᵀ, u = e_a + sign(β) e_b, ρ = |β|, β the coupling.
void	numerical::detail::solve_divide_leaf (const divide_node &node)
	Solves the leaf node by implicit QL: its eigenvalues ascending into its diagonal, their eigenvectors as the columns of its block of q.
void	numerical::detail::join_divide_node (const divide_node &node, merge_workspace &work)
	Joins the two solved halves of node by merge_rank_one, in work.
void	numerical::detail::unscale_divide_node (const divide_node &node)
	The eigenvalues of a solved node, scaled back, where its parent reads them.
void	numerical::detail::solve_divide_subtree (std::vector< divide_node > &nodes, std::size_t index, merge_workspace &work)
	Solves the block index of nodes and everything below it, each block right after its halves.
void	numerical::detail::divide_and_conquer (double *d, const double *e, std::size_t n, buffer< double > &q)
	The eigenvalues, ascending, of the symmetric tridiagonal with diagonal d[0…n-1] and subdiagonal e[0…n-2] into d, and its eigenvectors as the columns of the row-major q (n × n), by Cuppen's divide and conquer: the tridiagonal is two halves plus a rank-one coupling, each half is diagonalised the same way (by implicit QL up to divide_leaf), and merge_rank_one joins them, each half into its own diagonal block of q, so the joins copy no half and leave the blocks off the diagonal as they found them until the join of their two halves fills them.
template<typename T >
void	numerical::detail::multiply_real_left (const buffer< double > &left, buffer< T > &zt, std::size_t n)
	zt ← left · zt, n × n, row-major, for any other element type: the real factor converted.
void	numerical::detail::multiply_real_left (const buffer< double > &left, buffer< double > &zt, std::size_t n)
	zt ← left · zt, n × n, row-major, for a real left factor.
void	numerical::detail::multiply_real_left (const buffer< double > &left, buffer< std::complex< double > > &zt, std::size_t n)
	zt ← left · zt for a complex zt: the real and imaginary parts by a real product each.
template<typename T >
buffer< T >	numerical::detail::tridiagonal_divide_and_conquer (buffer< double > &d, buffer< double > &e, const tridiagonal_reduction< T > &reduction, std::size_t n)
	The eigenvalues of the tridiagonal (d, e) of reduction into d, and the matrix's eigenvectors as the rows of the result, by divide and conquer (LAPACK's dstedc, then dormtr): most of the work in matrix products. Up to divide_leaf, tridiagonal_ql on Mᵀ itself.
double	numerical::detail::jacobi_off_diagonal (const buffer< double > &a, std::size_t n)
	The Frobenius norm of the strict off-diagonal of the n × n a.
void	numerical::detail::jacobi_rotate (buffer< double > &a, buffer< double > &v, std::size_t n, std::size_t p, std::size_t q)
	One Jacobi rotation of the n × n a zeroing a[p][q], accumulated into v.
void	numerical::detail::cyclic_jacobi (buffer< double > &a, buffer< double > &v, std::size_t n, double frobenius, std::size_t max_sweeps)
	Cyclic Jacobi on the full symmetric a, the rotations accumulated into v.
std::vector< std::size_t >	numerical::detail::ascending_order (const buffer< double > &diagonal)
	The indices of diagonal by ascending value.
template<typename T >
buffer< T >	numerical::detail::columns_by_value (const buffer< T > &zt, const std::vector< std::size_t > &order, std::size_t n)
	The n × n row-major matrix whose column c is row order[c] of zt (the eigenvectors as rows): by tiles, the rows of each on its own core. Element by element, each write was a row apart, and the copy of order 1024 took a tenth of an eigh on six x86-64 cores.
template<typename T >
void	numerical::detail::one_sided_jacobi (buffer< T > &w, buffer< T > &v, std::size_t n, std::size_t max_sweeps)
	One-sided (Hestenes) Jacobi on the columns of the n × n row-major w: plane rotations, accumulated into v, until every pair of columns is orthogonal to rounding. The column norms are then the singular values of the input, each to high relative accuracy.
template<typename T >
void	numerical::detail::form_bidiagonal_bases (const std::vector< buffer< T > > &left, std::vector< buffer< T > > right, const buffer< T > &left_phase, const buffer< T > &right_phase, std::size_t n, buffer< T > &ut, buffer< T > &vt)
	The bases of a bidiagonal reduction from its reflectors: U_b = (H₀ H₁ …) diag(conj left phases) and V_b = (H'₀ H'₁ …) diag(right phases), each formed by blocks of reflectors through the matrix product kernel, as QR's Q is, or a reflector at a time below bidiagonal_basis_crossover.
template<typename T >
void	numerical::detail::bidiagonalize_unblocked (buffer< T > &b, std::size_t n, buffer< double > &d, buffer< double > &e, buffer< T > &ut, buffer< T > &vt, bool form)
	Reduces the square n × n row-major b to real upper-bidiagonal form, B = U_b D V_bᴴ with D real, by Householder reflections alternately from the left and the right (Golub–Kahan).
template<typename T >
T	numerical::detail::householder_in_place (buffer< T > &x)
	The unit Householder vector v with (I − 2 v vᵀ) x = β e₁, in place of x, for a real x.
template<typename T > requires std::is_floating_point_v<T>
void	numerical::detail::bidiagonalize (buffer< T > &b, std::size_t n, buffer< double > &d, buffer< double > &e, buffer< T > &ut, buffer< T > &vt, bool form=true)
	Reduces the real square n × n row-major b to upper-bidiagonal form, B = U_b D V_bᵀ, by panels of bidiagonal_block reflections from each side, as LAPACK's dgebrd and dlabrd: the matrix a panel's later steps see is B − V Yᵀ − X Uᵀ (V, U the reflectors, Y, X what they move), applied to the one column and the one row each step needs; the trailing block takes it after the panel, through the matrix product kernel. The bidiagonal keeps its signs, which the iterations that follow take.
template<typename T > requires (!std::is_floating_point_v<T>)
void	numerical::detail::bidiagonalize (buffer< T > &b, std::size_t n, buffer< double > &d, buffer< double > &e, buffer< T > &ut, buffer< T > &vt, bool form=true)
	bidiagonalize_unblocked, for a complex or any other T the panels do not take.
template<typename T >
void	numerical::detail::bidiagonal_qr (buffer< double > &d, buffer< double > &e, buffer< T > &ut, buffer< T > &vt, std::size_t n, std::size_t max_steps)
	The singular values of the real upper bidiagonal (d, e) by implicit-shift QR with deflation (LINPACK's dsvdc, the Golub–Kahan–Reinsch iteration), the rotations accumulated into the rows of ut and vt.
template<typename T >
matrix< T >	numerical::detail::adjoint (const matrix< T > &a)
	The conjugate transpose of a.
double	numerical::detail::singular_secular_root (const buffer< double > &d, const buffer< double > &z, std::size_t i, std::size_t max_iterations, double *below, double *above)
	Singular value i of the arrow M whose first row is z and whose diagonal is d (d₀ = 0, d strictly ascending): root σ of 1 + Σⱼ zⱼ²/(dⱼ² − σ²) = 0, as LAPACK's dlasd4.
void	numerical::detail::rotate_columns (buffer< double > &m, std::size_t cols, std::size_t a, std::size_t b, double c, double s)
	Columns a and b of the row-major m (cols per row) rotated: a ← c a + s b, b ← c b − s a.
arrow_workspace &	numerical::detail::thread_arrow_workspace ()
	The calling thread's arrow_workspace, kept from one decomposition to the next. A join uses it whole on one thread and nests no other join; release_large_arrow_workspace lets a large buffer go after.
void	numerical::detail::release_large_arrow_workspace (arrow_workspace &work)
	Releases the buffers of work longer than kept_workspace.
void	numerical::detail::merge_arrow (buffer< double > &d, buffer< double > z, double alpha, double beta, buffer< double > &u, buffer< double > &v, std::size_t n, std::size_t m, std::size_t split, arrow_workspace &work)
	The singular value decomposition of the arrow M (first row z, diagonal d, d₀ = 0), carried into the bases u (n × n) and v (m × m, its first n columns the arrow's; a last, null column when m = n + 1), as LAPACK's dlasd1–dlasd3: on return d holds the singular values ascending and u, v the singular vectors, v's null column last.
std::size_t	numerical::detail::plan_bidiagonal (std::vector< bidiagonal_node > &nodes, const double *d, const double *e, std::size_t n, std::size_t extra, std::size_t depth)
	Appends the block (d, e, n, extra) and everything below it to nodes, parents before children, each scaled to work near unit size by a power of two (exactly): the secular equation works on squares, which overflow past 1e154.
void	numerical::detail::start_bidiagonal_node (bidiagonal_node &node)
	Its singular values and vectors into node's, from scratch: u zero, v the identity.
void	numerical::detail::solve_bidiagonal_leaf (bidiagonal_node &node)
	Decomposes the leaf node by implicit-shift QR.
void	numerical::detail::join_bidiagonal_node (bidiagonal_node &node, const bidiagonal_node &one, const bidiagonal_node &two, arrow_workspace &work)
	Joins the two decomposed halves of node, one and two, by merge_arrow, node started (start_bidiagonal_node).
void	numerical::detail::unscale_bidiagonal_node (bidiagonal_node &node)
	The singular values of a solved node scaled back.
void	numerical::detail::join_bidiagonal_halves (std::vector< bidiagonal_node > &nodes, std::size_t index, arrow_workspace &work)
	Joins the halves of nodes' block index in work, then lets them go.
void	numerical::detail::solve_bidiagonal_subtree (std::vector< bidiagonal_node > &nodes, std::size_t index, arrow_workspace &work)
	Decomposes the block index of nodes and everything below it, each block right after its halves.
void	numerical::detail::bidiagonal_divide (const double *d, const double *e, std::size_t n, std::size_t extra, buffer< double > &sigma, buffer< double > &u, buffer< double > &v)
	The singular value decomposition of the upper bidiagonal n × (n + extra) (extra 0 or 1) with diagonal d[0…n-1] and superdiagonal e[0…n+extra-2] (e[n-1] the last column's entry when extra is 1), by divide and conquer, as LAPACK's dbdsdc: sigma ascending, u (n × n) and v (m × m, m = n + extra, the null vector last) row-major, their columns the singular vectors.
template<typename T >
void	numerical::detail::bidiagonal_divide_and_conquer (buffer< double > &d, buffer< double > &e, buffer< T > &ut, buffer< T > &vt, std::size_t n)
	The singular values of the upper bidiagonal (d, e) into d, descending, the singular vectors accumulated into the rows of ut and vt as bidiagonal_qr accumulates them, by divide and conquer (bidiagonal_divide): most of the work in matrix products.
std::size_t	numerical::detail::numerical_rank (const vector< double > &sigma, std::size_t rows, std::size_t cols, rank_cutoff cutoff={})
	How many of the singular values sigma of an rows × cols matrix exceed cutoff: its numerical rank (the first that many, where sigma is descending).
template<typename T >
T	numerical::detail::scaled_by_two (const T &x, int exponent)
	x times two to the exponent, exactly where the result is a normal number (no power of two is formed, so none overflows).
template<typename U >
std::complex< U >	numerical::detail::scaled_by_two (const std::complex< U > &x, int exponent)
	scaled_by_two for a complex number, part by part.
template<typename T >
buffer< T >	numerical::detail::upper_inverse (const buffer< T > &r, std::size_t n)
	R⁻¹, upper triangular, of the n × n upper-triangular r (row-major): back substitution against the identity in parallel chunks of columns, each row's terms in registers (subtract_products). A chunk's rows past its last column stay zero; a zero on the diagonal leaves infinities or NaN.
template<typename T >
std::optional< matrix< T > >	numerical::detail::lstsq_by_qr (const matrix< T > &a, const matrix< T > &b, rank_cutoff cutoff, std::pair< buffer< T >, int > *const r=nullptr)
	The least-squares solution of A X ≈ B by a QR factorisation, of A, or of Aᴴ for a wide A (the minimum-norm solution), where it is the one the singular value decomposition gives: with A divided by a power of two near its largest magnitude, the bound 1 / ‖R⁻¹‖_F under its smallest singular value must stay qr_rank_margin times above cutoff of the bound ‖R‖_F over its largest, so that every singular value is far above the cutoff and the rank full (a zero on R's diagonal makes the first bound zero, or NaN, which fails). Else nothing. A fraction of the decomposition's work, the rest in matrix products.
template<typename T >
buffer< double >	numerical::detail::square_singular_values (buffer< T > b, std::size_t n)
	The singular values of the square n × n row-major b, descending, without its vectors: its Golub–Kahan bidiagonal and the implicit-shift QR iteration on it, no basis formed nor rotated.
template<typename T >
matrix< T >	numerical::detail::lstsq_from_svd (const svd_decomposition< T > &s, std::size_t rank, const matrix< T > &b)
	The minimum-norm least-squares solution of A X ≈ B from the singular value decomposition s of A (m × n), its first rank triplets: X = V Σ⁺ Uᴴ B.
template<typename T >
vector< T >	numerical::lstsq (const matrix< T > &a, const vector< T > &b, rank_cutoff cutoff={})
	The minimum-norm least-squares solution of A x ≈ b, as NumPy's lstsq gives it.
template<typename T >
matrix< T >	numerical::lstsq (const matrix< T > &a, const matrix< T > &b, rank_cutoff cutoff={})
	The minimum-norm least-squares solution of A X ≈ B for every column of B at once, as numpy.linalg.lstsq with a matrix right-hand side: one SVD of A, X = V Σ⁺ Uᴴ B; or, where a bound on the condition number shows the rank full with room to spare, the same solution to rounding by one QR factorisation (detail::lstsq_by_qr).
template<typename T >
vector< double >	numerical::singular_values (const matrix< T > &a)
	The singular values of a, descending, without its singular vectors, as numpy.linalg.svd(a, compute_uv=False): the Householder R of a tall A (of Aᴴ when wide), its Golub–Kahan bidiagonal and the implicit-shift QR iteration on it, no basis formed nor rotated. Each value to an absolute accuracy of the largest times the machine epsilon, as svd_decomposition's.
template<typename T >
vector< double >	numerical::symmetric_eigenvalues (const matrix< T > &a)
	The eigenvalues of the real symmetric or complex Hermitian a, ascending, without its eigenvectors, as numpy.linalg.eigvalsh: its Householder tridiagonal (no basis formed) and the implicit QL iteration on it, no rotation kept. Only the lower triangle is read, and of the diagonal its real part. Each value to an absolute accuracy of the largest magnitude times the machine epsilon, as symmetric_eigen's.
template<typename T >
least_squares_solution< T >	numerical::least_squares (const matrix< T > &a, const matrix< T > &b, rank_cutoff cutoff={})
	lstsq, with the residuals, the rank and the singular values numpy.linalg.lstsq returns beside the solution: the residual sums of squares are given only where the solution is unique and need not be exact, A tall (m > n) and of full column rank, as NumPy gives them.
template<typename T >
matrix< T >	numerical::detail::hermitian_from_lower (const matrix< T > &a)
	The real symmetric or complex Hermitian matrix whose lower triangle is the square a's, and of its diagonal the real part.
template<typename T >
matrix< T >	numerical::detail::inverse_on_kept (const vector< double > &sigma, const matrix< T > &right, const matrix< T > &left, rank_cutoff cutoff)
	Σ_d right_d (1 / σ_d) left_dᴴ over the columns d of right (cols × k) and left (rows × k) whose sigma is above cutoff (of the largest magnitude in sigma): a pseudo-inverse from a singular value decomposition (right = V, left = U) or from an eigendecomposition (V both, σ the eigenvalues, signed). Through the matrix product kernel: summed term by term it ran five times as long as NumPy's at order 512.
template<typename T >
vector< double >	numerical::detail::hermitian_singular_values (const matrix< T > &a)
	The singular values of the real symmetric or complex Hermitian matrix whose lower triangle is a's: of the matrix made whole, for a scalar the eigendecompositions do not take (they compute in double).
template<typename T >
vector< double >	numerical::detail::eigenvalue_magnitudes (const matrix< T > &a)
	The magnitudes of the eigenvalues of the real symmetric or complex Hermitian a, its lower triangle read: its singular values, unordered.
vector< double >	numerical::detail::hermitian_singular_values (const matrix< double > &a)
	hermitian_singular_values in double: the magnitudes of the eigenvalues, cheaper.
vector< double >	numerical::detail::hermitian_singular_values (const matrix< std::complex< double > > &a)
	hermitian_singular_values in std::complex<double>: the magnitudes of the eigenvalues, cheaper.
template<typename T >
matrix< T >	numerical::detail::hermitian_pseudo_inverse (const matrix< T > &a, rank_cutoff cutoff)
	The pseudo-inverse of the real symmetric or complex Hermitian matrix whose lower triangle is a's: from the singular value decomposition of the matrix made whole, for a scalar the eigendecompositions do not take.
matrix< double >	numerical::detail::hermitian_pseudo_inverse (const matrix< double > &a, rank_cutoff cutoff)
	hermitian_pseudo_inverse in double: V Λ⁺ Vᵀ from the eigendecomposition, cheaper.
matrix< std::complex< double > >	numerical::detail::hermitian_pseudo_inverse (const matrix< std::complex< double > > &a, rank_cutoff cutoff)
	hermitian_pseudo_inverse in std::complex<double>: V Λ⁺ Vᴴ from the eigendecomposition, cheaper.
template<typename T >
std::size_t	numerical::matrix_rank (const matrix< T > &a, rank_cutoff cutoff={}, bool hermitian=false)
	The numerical rank of a, as numpy.linalg.matrix_rank counts it: its singular values above cutoff.
template<typename T >
matrix< T >	numerical::pinv (const matrix< T > &a, rank_cutoff cutoff={}, bool hermitian=false)
	The Moore–Penrose pseudo-inverse of the m × n a (n × m), as numpy.linalg.pinv gives it: V Σ⁺ Uᴴ, a singular value at or under cutoff counting as zero.

Variables
constexpr std::size_t	numerical::detail::lu_block {64}
	The columns of one panel of the blocked LU: the trailing update, a product of this depth, is deep enough for the packed kernel while the unblocked panel stays a small part of the work.
constexpr std::size_t	numerical::detail::lu_small_order {256}
	The order up to which the LU takes panels of lu_small_block columns: there the panels, a column at a time on one thread, are most of the time, and narrower ones halve their work (arm64, order 300 on every core: 14 % slower at 32 columns than at 64).
constexpr std::size_t	numerical::detail::lu_small_block {32}
	The columns of one panel of the LU of an order up to lu_small_order.
constexpr std::size_t	numerical::detail::lu_lookahead_columns {1024}
	The trailing columns at most past which the LU factorises its next panel beside the update of the rest (a look-ahead): a wider update runs well as one product, and cut in blocks beside the panel it ran longer (order 2048 on six x86-64 cores, 9 %); a narrower one leaves threads idle while the panel runs.
constexpr std::size_t	numerical::detail::lu_update_speedup {gemm_kernel_speedup}
	The multiply-adds of an LU's trailing update one of the tile's dependent ones stands for, as threads are weighed (gemm_kernel_speedup): the AMX tile's, on arm64; the AVX tiles' quarter of it on x86-64, where the update of order 256 waited on one thread (six cores: 690 → 622 µs, 1024: 12.5 → 12.1 ms).
constexpr std::size_t	numerical::detail::narrow_solve {16}
	The most right-hand sides an LU solve takes row by row rather than by blocks through the product kernel.
constexpr std::size_t	numerical::detail::slice_columns {64}
	The fewest columns an LU solve gives a thread of its own: a narrower slice starves the product kernel (arm64, order 256 and 64 right-hand sides in eight slices: 112 → 194 µs; 48 columns, order 300 and 300 right-hand sides: 8 % slower than 64).
constexpr std::size_t	numerical::detail::solve_leaf {16}
	The columns of one leaf of the recursive LU panel.
constexpr std::size_t	numerical::detail::inverse_band {16}
	The rows an LU inverse writes its identity into, and permutes the columns of, on one thread at a time.
constexpr std::size_t	numerical::detail::memory_element_weight {4}
	The weight of an element a pass bound by the memory writes or moves, against a multiply-add, as the BLAS routines' transposing copies (weighed as one, an LU inverse of order 300 permuted its columns on one thread).
constexpr std::size_t	numerical::detail::svd_band {32}
	The rows of an SVD's input one thread scales and copies at a time.
constexpr std::size_t	numerical::detail::zeroed_rows {32}
	The rows of a Cholesky factor whose upper part one thread zeroes at a time.
constexpr std::size_t	numerical::detail::cholesky_block {64}
	The columns of one panel of the blocked Cholesky factorisation (as lu_block).
constexpr std::size_t	numerical::detail::cholesky_diagonal_rows {4}
	The rows of a Cholesky diagonal block factorised together, one chain of sums each (one x86-64 AVX2 core, order 300: the diagonal blocks 128 → 65 µs; eight rows did no better).
constexpr std::size_t	numerical::detail::cholesky_together_width {16}
	The narrowest Cholesky diagonal block whose rows go together: in a block of 8, the groups' loops cost more than the chains save (order 8, one arm64 core: 0.21 → 0.23 µs).
constexpr std::size_t	numerical::detail::cholesky_one_panel {48}
	The order below which the Cholesky factorisation is one panel: no trailing product, whose packing the few panels of a small order do not repay (order 17 on one thread: 0.8 µs whole, 1.3 µs in panels of 16). Panels of 16 win from order 40 on x86-64 and from 48 on arm64.
constexpr std::size_t	numerical::detail::cholesky_rows {64}
	The rows of one chunk of a Cholesky panel's rows below it, solved together along each column.
constexpr std::size_t	numerical::detail::cholesky_update_columns {64}
	The fewest columns of the trailing block one Cholesky update's product takes: each product packs every row below its columns again, a panel's width apart the packing outweighed the products on the AMX (M1 Pro, order 128: 61 → 53 µs), and the triangle a wider block computes past the diagonal is wasted.
constexpr std::size_t	numerical::detail::cholesky_update_rows {256}
	The rows of one tile of a Cholesky update: the tiles of a panel's update, each one product on one thread, share the threads in one parallel loop, where a parallel product for each block of columns started and joined the threads some thirty times a panel.
template<typename T >
constexpr double	numerical::detail::rescue_below {std::is_same_v<value_range_t<T>, float> ? 0x1p-60 : 0x1p-500}
	The magnitude under which rescue_factor scales a T up: about the square root of the smallest normal number of its range.
constexpr double	numerical::detail::plain_squares_low {0x1p-900}
	The smallest plain sum of squares plain_or_at_scale keeps: a term that underflowed is below 2⁻¹⁰²², so a sum this large loses at most count · 2⁻¹²² of itself to them.
constexpr double	numerical::detail::plain_squares_high {0x1p900}
	The largest plain sum of squares plain_or_at_scale keeps: no term of it overflowed.
constexpr std::size_t	numerical::detail::finite_chunk {std::size_t {1} << 15U}
	The elements all_finite_shared gives a thread at a time.
constexpr std::size_t	numerical::detail::qr_block {32}
	The reflectors of one panel of the blocked QR.
constexpr std::size_t	numerical::detail::qr_transpose_tile {32}
	The side of the tiles in which the QR transposes (transpose_into): a tile's rows and columns both stay in the cache, where a whole column read across rows misses on every element.
constexpr std::size_t	numerical::detail::jacobi_sweeps {60}
	The sweeps a Jacobi method may take: it converges quadratically, in under 15 on every matrix measured; the bound turns a non-convergence into an error instead of a hang.
constexpr std::size_t	numerical::detail::corrected_product_least {512}
	The rows from which corrected_hermitian_product shares its coefficients and bands with the threads: on six x86-64 AVX2 cores from 256 (below, eigh of order 128 took 13 % longer), on an M1 Pro from 512 (from 256, eigh of order 512 on eight cores took 3 to 9 % longer).
constexpr std::size_t	numerical::detail::tridiagonal_block {32}
	The reflections of one panel of the blocked tridiagonal reduction.
constexpr std::size_t	numerical::detail::tridiagonal_band {512}
	The rows of each band of the tridiagonal reduction's trailing update: a band reaches its diagonal, so a wider one computes more of the upper triangle the reduction never reads, a narrower one runs the product kernel on smaller products.
constexpr std::size_t	numerical::detail::ql_steps_per_value {30}
	The implicit QL steps the tridiagonal eigensolver may take per eigenvalue, as LAPACK's dsteqr bounds them (30 n in all): it converges cubically, in two or three per eigenvalue.
constexpr std::size_t	numerical::detail::divide_leaf {32}
	The order up to which tridiagonal_divide_and_conquer diagonalises by implicit QL.
constexpr std::size_t	numerical::detail::secular_iterations {100}
	The iterations one root of the secular equation may take: it converges in a few; the bound turns a non-convergence into an error instead of a hang.
constexpr std::size_t	numerical::detail::lowner_block {32}
	The poles a task of merge_rank_one's Löwner products takes, and the side of its tiles of u.
constexpr std::size_t	numerical::detail::first_touch_chunk {std::size_t {1} << 15}
	The elements a task of workspace_of touches first: 64 pages of 4 KiB, so a buffer of a few megabytes spreads over every core.
constexpr std::size_t	numerical::detail::kept_workspace {std::size_t {1} << 18}
	The most elements a buffer of a thread's merge_workspace keeps from one decomposition to the next: 2 MiB, every buffer an eigh of order 512 needs, so a thread keeps at most 8 MiB.
constexpr double	numerical::detail::jacobi_tiny {std::numeric_limits<double>::min() / std::numeric_limits<double>::epsilon()}
	The sum of squares under which one-sided Jacobi judges a pair of columns over their norms: below it a square has lost digits to underflow (DBL_MIN / ε).
constexpr double	numerical::detail::jacobi_large {1e150}
	The |ζ| past which one-sided Jacobi takes the rotation's tangent as 1 / (2|ζ|): ζ² stays finite under it.
constexpr std::size_t	numerical::detail::svd_qr_steps_per_value {30}
	The QR steps the bidiagonal singular value iteration may take per value: it converges in two or three; the bound turns a non-convergence into an error instead of a hang.
constexpr std::size_t	numerical::detail::bidiagonal_crossover {128}
	The order below which a real bidiagonal reduction goes a reflection at a time, as LAPACK's dgebrd falls back to dgebd2: the panels' products are too small to pay for packing. Order 64: 33 % shorter on six x86-64 cores, 7 % on arm64; order 160: 30 % longer on arm64, where x86-64 still gains 8 %.
constexpr std::size_t	numerical::detail::bidiagonal_basis_crossover {64}
	The order below which the bases of a bidiagonal reduction are formed a reflector at a time: the blocks' packing and triangular factor cost more than they save. Order 32: 11 % shorter on six x86-64 cores, 12 % on arm64; order 96: 7 % longer on both.
constexpr std::size_t	numerical::detail::bidiagonal_block {32}
	The rows from which a step of the blocked bidiagonal reduction shares its coefficients and corrections with the threads: on x86-64 (six AVX2 cores, the reduction of order 512 −5 to −14 %, 1024 −12 %); never on arm64, where its serial part is smaller and the shared step gained nothing (M1 Pro, eight cores, order 1024 0.97 to 1.05 times as long).
constexpr std::size_t	numerical::detail::bidiagonal_product_weight {1}
	The weight of a multiply-add of the bidiagonal reduction's matrix-vector products, which read the trailing block once for each column: on x86-64 a thread shares them from a quarter of the work (six cores, singular values of order 256: 8.1 → 7.3 ms, 512: 35 → 30); on arm64 that slowed order 256 by 10 %.
constexpr std::size_t	numerical::detail::bidiagonal_leaf {32}
	The order up to which bidiagonal_divide decomposes by implicit-shift QR.
constexpr double	numerical::detail::qr_rank_margin {100.0}
	How far above the rank cutoff a bound on the smallest singular value must stay for lstsq_by_qr to answer: nearer, the singular value decomposition decides the rank.

Detailed Description

Matrix factorisations and the algorithms built on them.

numerical::lu_decomposition factorises a square matrix once (Gaussian elimination with partial pivoting — an in-place P A = L U), then answers any number of queries cheaply: numerical::lu_decomposition::solve a system, its numerical::lu_decomposition::determinant, or its numerical::lu_decomposition::inverse. Free numerical::solve, numerical::determinant and numerical::inverse are one-shot conveniences over it.

The scalar T must be a field (support division) — double, std::complex<double>, … — and the partial pivoting is numerically stable for well-conditioned systems. A singular matrix has a zero determinant; asking it to solve or invert throws numerical::singular_matrix_error.

Every factorisation here refuses an infinite or NaN entry where it starts, with a numerical::non_finite_error, since none of them means anything on one and the iterative ones would stop at once and return plausible values.

Definition in file decomposition.hpp.