|
NumericAL
A header-only C++20 lexer built on REAL
|
Classes | |
| struct | a_rows |
A as gemm_packed reads it by default: row-major, m × k, its rows lda apart, packed by the tile. More... | |
| struct | a_transposed |
A given as Aᵀ: row-major, k × m, its rows lda apart (C = Aᵀᵀ·B, no copy of the transpose). More... | |
| struct | 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 | axis_walk |
The extents and strides that walk a tensor of shape with axis axis taken out: outer positions before it, inner after it, the axis itself extent long at stride inner. More... | |
| struct | b_rows |
B as gemm_packed reads it by default: row-major, k × n, its rows ldb apart. More... | |
| struct | b_transposed |
B given as Bᵀ: row-major, n × k, its rows ldb apart (C = A·Bᵀᵀ, no copy of the transpose). More... | |
| struct | balancing |
How a matrix was balanced: swaps, in the order they were applied, then the scaling of the active block low … high − 1 by scale. More... | |
| struct | 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... | |
| struct | complex_operand |
A complex vector that rows are summed against, kept as (re, −im) and (im, re) pairs: a complex row, read as its parts in order ([complex.numbers]), against each is the real or the imaginary part of the sum, one dot_product each. More... | |
| struct | contraction_step |
One step of a contraction path: the operands at first and second of the current list, whose result goes to the list's end. More... | |
| struct | default_init_allocator |
The allocator of numerical::buffer. An element constructed with no value is default-initialised, so a double keeps whatever its memory held, where std::allocator zeroes it; one constructed with a value is constructed from it. A result written in full is then written once, by the threads that compute it, not zeroed first on one of them (a product of order 512 on eight arm64 threads spent a fifth of its time in that zeroing). More... | |
| struct | deflation_count |
| What early_deflation found: the eigenvalues left as shifts, and those deflated. More... | |
| struct | 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... | |
| struct | hessenberg_panel |
One panel of hessenberg_blocked, Q_panel = I − V T Vᴴ. More... | |
| struct | interleaved_components |
A complex matrix's eigenvector during its back-substitution, kept as (re, −im) and (im, re) pairs: a complex row, read as its parts in order ([complex.numbers]), against each is the real or the imaginary part of their product, one dot_product each. More... | |
| struct | 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 | one_based |
A row-major matrix read through indices from 1, (1, 1) its first entry. More... | |
| struct | packed_left |
A left operand packed once for several products by it (an LU's −L₂₁ against blocks of its U₁₂): the tile chosen and A packed for it (gemm_packed's layout of one panel), or A to read in place (plain). More... | |
| class | packing_storage |
| The storage of a product's packed operands, kept by the calling thread from one product to the next. More... | |
| struct | plane_rotation |
A plane rotation [c s; −s̄ c] and what it leaves of its first argument. More... | |
| struct | pool_job |
| One parallel_for handed to the pool: its chunks, the next one to take, and how many are done. Each run has its own, so a worker still leaving one run never touches the next. More... | |
| struct | qr_reflections |
The reflections of a Householder QR, Q left as its blocks of reflectors. More... | |
| struct | real_components |
| A real eigenvector during its back-substitution, the sums of rows against it by dot_product. More... | |
| struct | real_eigenvectors |
A real matrix's eigenvectors before their back-transformation: Xᵀ's real and imaginary parts. More... | |
| struct | real_operand |
| A real vector that rows are summed against by dot_product. More... | |
| struct | reflector_block |
The reflectors of one panel of the blocked QR, I − V T Vᴴ as LAPACK's dlarft forms it. More... | |
| struct | 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 | row_operation |
| One operation of tridiagonal_ql or bidiagonal_qr on the rows of a vector matrix, recorded to be applied later. More... | |
| struct | scaled_quotient |
A quotient of divided and the factor its dividend was scaled by: a type of its own, where a std::pair of them made GCC on arm64 print, at every caller, a note on its ABI change of GCC 10.1. More... | |
| struct | schur_target |
Where a QR sweep applies its reflections: the rows and columns it updates, and Z. More... | |
| struct | secular_solution |
| A root of a secular equation: the pole it was sought from, and its offset from that pole. More... | |
| struct | split_components |
| A real matrix's complex eigenvector during its back-substitution, its real and imaginary parts apart: a row of the real Schur form against it is two dot_product. More... | |
| struct | split_operands |
The real operands of gemm_split, kept by each thread from one product to the next (split_storage): mapping their pages anew took a third to a half of a complex product of order 64 to 256 on one x86-64 core (six AVX2 cores, all threads: order 128 0.52 → 0.17 ms), 15–20 % of one of order 256 on arm64. More... | |
| struct | standard_block |
| The rotation standardize_block found, and the eigenvalues of the block it standardised. More... | |
| struct | thin_qr |
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... | |
| class | thread_pool |
| The threads every kernel shares, started once: a parallel_for costs a wake-up, not a thread start. More... | |
| struct | tridiagonal_reduction |
The unitary M = Q Dᴴ of a reduction to real tridiagonal form, A = M T Mᴴ, kept as its factors. More... | |
| struct | vector_tile_kernel |
The vector path's tile: gemm_mr × gemm_nr elements of C in registers (gemm_kernel), each product rounded before its sum unless the compiler contracts c + a·b into a fused multiply-add. More... | |
| struct | worker_bell |
| A worker's own wake-up: its core, whether it sleeps, and what wakes it. More... | |
| struct | workspace_buffer |
One buffer of merge_workspace: its storage and how many elements it holds. More... | |
Concepts | |
| concept | split_complex |
A complex number of float or double, which gemm_split multiplies through the real kernel and the kernels take apart into its real and imaginary parts. | |
Typedefs | |
| template<typename T > | |
| using | 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 | 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 | 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 | 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. | |
| using | lapack_index = std::ptrdiff_t |
| A signed index from 1, as LAPACK's: the multishift sweep starts its bulges in phantom columns above the matrix, before its first row. | |
| template<typename T > | |
| using | gemm_tile = std::array< T, gemm_mr *gemm_nr< T > > |
| template<term_rounding Mode> | |
| using | rounding_is = std::integral_constant< term_rounding, Mode > |
The rounding mode Mode as a type, which picks rounded_term's overload. | |
| using | chunk_runner = void(*)(const void *body, std::size_t begin, std::size_t end) |
A chunk of a parallel_for, type-erased: runs body(i) for i in [begin, end). | |
| using | side_runner = void(*)(const void *side) |
| A task the caller of a run does before it takes chunks (parallel_for_beside), type-erased. | |
Enumerations | |
| enum class | term_rounding : unsigned char { compiler , fused , x86_fused } |
How numerical::detail::gemm_rows rounds each term sum + a·b, as the tile of the path in use does: as the compiler rounds it (the vector tile's way), one fused multiply-add (arm64's NEON and AMX tiles), or one fused multiply-add of AVX2's unit (the x86-64 tiles), in code built for any x86-64 processor. More... | |
| enum class | left_kernel : unsigned char { plain , vector , neon , amx , avx512 , avx2 } |
The tile a product's left operand is packed for (numerical::detail::pack_left): plain where it is not packed. More... | |
Functions | |
| bool | amx_usable () noexcept |
| Whether this chip's AMX can be used: known family and a passing self-test, decided once. | |
| template<typename T > | |
| bool | operator== (const std::vector< T, default_init_allocator< T > > &buffer, const std::vector< T > &other) |
| Element-wise equality of a numerical::buffer and a std::vector, found by argument-dependent lookup through the allocator: code that compared a tensor's elements with a std::vector keeps doing so. | |
| template<typename T > | |
| T | 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 > | pivot_divisor (const std::complex< R > &pivot) |
| pivot_divisor of a complex pivot: its reciprocal. | |
| template<typename T > | |
| T | 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 > | 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 | 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 | 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 | 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 | 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 | conj_of (const T &x) |
| The conjugate of a real (or general) scalar — the value itself. | |
| template<typename U > | |
| std::complex< U > | conj_of (const std::complex< U > &x) |
| The conjugate of a complex scalar. | |
| template<typename T > | |
| T | 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 | real_part (const T &x) |
The real part of a real (or general) scalar as a double. | |
| template<typename U > | |
| double | real_part (const std::complex< U > &x) |
| The real part of a complex scalar. | |
| template<typename T > | |
| double | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | scaled_norm (const buffer< T > &x) |
| scaled_norm of a whole vector. | |
| double | 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 | is_finite_entry (const T &x) |
Whether x is finite. | |
| template<typename U > | |
| bool | is_finite_entry (const std::complex< U > &x) |
Whether x is finite: both of its parts. | |
| template<typename T > | |
| void | require_finite (const matrix< T > &a, const char *what) |
| The one non-finite policy, checked where every decomposition starts. | |
| bool | 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 > | 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 | 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 | require_finite (const matrix< double > &a, const char *what) |
require_finite for double, its elements contiguous: all_finite. | |
| void | 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 | 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 | 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 | 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 | 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 | 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 | 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 > | |
| void | 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 > | 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 | 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 | 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 | 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 | 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 | 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 | 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 > > | 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 | 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 | 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 > > | 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 > | 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 | 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 > | 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 | 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 > requires std::is_floating_point_v<T> | |
| void | 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 | 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 | 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 > | 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 | 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 | 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 | 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 | 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 | 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 | 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 | real_subdiagonal (double t, double, double &) |
| A real subdiagonal entry of the tridiagonal: kept with its sign (the rotations take any sign). | |
| double | 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 > | 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 > | 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 | 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 > | 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 | 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 | 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 | 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 | 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 | 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 * | 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 & | 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 | release_large_workspace (merge_workspace &work) |
Releases the buffers of work longer than kept_workspace. | |
| void | 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 | 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 | 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 | join_divide_node (const divide_node &node, merge_workspace &work) |
Joins the two solved halves of node by merge_rank_one, in work. | |
| void | unscale_divide_node (const divide_node &node) |
The eigenvalues of a solved node, scaled back, where its parent reads them. | |
| void | 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 | 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 | 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 | 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 | 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 > | 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 | 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 | 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 | 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 > | ascending_order (const buffer< double > &diagonal) |
The indices of diagonal by ascending value. | |
| template<typename T > | |
| buffer< T > | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 > | adjoint (const matrix< T > &a) |
The conjugate transpose of a. | |
| double | 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 | 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 & | 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 | release_large_arrow_workspace (arrow_workspace &work) |
Releases the buffers of work longer than kept_workspace. | |
| void | 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 | 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 | start_bidiagonal_node (bidiagonal_node &node) |
Its singular values and vectors into node's, from scratch: u zero, v the identity. | |
| void | solve_bidiagonal_leaf (bidiagonal_node &node) |
Decomposes the leaf node by implicit-shift QR. | |
| void | 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 | unscale_bidiagonal_node (bidiagonal_node &node) |
The singular values of a solved node scaled back. | |
| void | 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 | 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 | 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 | 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_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 | 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 > | 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 > | 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 > > | 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 > | 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 > | 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 > | |
| matrix< T > | 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 > | 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 > | 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 > | 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 > | hermitian_singular_values (const matrix< double > &a) |
hermitian_singular_values in double: the magnitudes of the eigenvalues, cheaper. | |
| vector< double > | 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 > | 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 > | 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 > > | 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. | |
| double | modulus_1 (double x) |
|x| of a real x: what modulus_1 of a complex one stands for. | |
| double | modulus_1 (const std::complex< double > &x) |
|re x| + |im x|, the cheap modulus LAPACK's complex routines judge with. | |
| template<typename A , typename B > | |
| auto | times (const A &a, const B &b) |
a · b. | |
| std::complex< double > | times (const std::complex< double > &a, const std::complex< double > &b) |
a · b of two complex numbers by the schoolbook formula, without C Annex G's recovery of infinities, which a large translation unit leaves out of line (a call per product: complex eig of order 64 took 2.7 times as long in the Python binding as alone). The entries here are finite. | |
| double | largest_part (double x) |
|x| of a real x. | |
| double | largest_part (const std::complex< double > &x) |
max(|re x|, |im x|): within √2 of |x|, without its square root. | |
| template<typename T > | |
| void | exchange (buffer< T > &a, std::size_t n, std::size_t i, std::size_t j) |
Exchanges rows i and j, then columns i and j, of the row-major n × n a: a similarity. | |
| template<typename T > | |
| balancing | balance (buffer< T > &a, std::size_t n) |
Balances the row-major n × n a in place, as LAPACK's ‘dgebal('B’)`: rows with no off-diagonal entry in the active block go to its bottom and columns with none to its top, each isolating an eigenvalue on the diagonal; then each row and column of the remaining block is scaled by a power of two until their 2-norms are within a factor of 2, which leaves every entry exact. | |
| template<typename T > | |
| buffer< T > | hessenberg_unblocked (buffer< T > &a, std::size_t n, std::size_t low, std::size_t high, bool form) |
Reduces the active block low … high − 1 of the balanced row-major n × n a to upper Hessenberg form by Householder reflections one column at a time, as LAPACK's dgehd2: A = Q H Qᴴ, the entries below the subdiagonal set to zero. | |
| bool | already_reduced (double, double rest) |
Whether (α, x), of tail norm rest, is already (β, 0) with β real: for a real one, its tail zero. | |
| bool | already_reduced (const std::complex< double > &alpha, double rest) |
The same for a complex one: its tail zero and α real. | |
| template<typename T > | |
| T | reflector_tau (T &alpha, T *x, std::size_t count) |
The Householder reflection H = I − τ v vᴴ, v = (1, v₁…), with Hᴴ (α, x) = (β, 0) and β real, as LAPACK's zlarfg (dlarfg for a real T): alpha becomes β and the count entries at x become v₁…. | |
| real_operand | operand_of (const double *v, std::size_t stride, std::size_t count) |
The count entries stride apart at v as a real_operand. | |
| complex_operand | operand_of (const std::complex< double > *v, std::size_t stride, std::size_t count) |
| The same as a complex_operand. | |
| template<typename T > | |
| buffer< T > | hessenberg_blocked (buffer< T > &a, std::size_t n, std::size_t low, std::size_t high, bool form) |
hessenberg_unblocked by panels of hessenberg_block columns, as LAPACK's dgehrd and dlahr2: a panel's reflectors are formed one column at a time, each column first updated by those before it through Y = A V T; the rest of the matrix then takes the panel from both sides through the matrix product kernel, A ← (I − V Tᴴ Vᴴ)(A − Y Vᴴ). Its arguments and answer are hessenberg_unblocked's. | |
| template<typename T > | |
| buffer< T > | hessenberg (buffer< T > &a, std::size_t n, std::size_t low, std::size_t high, bool form) |
| hessenberg_blocked from an active order of hessenberg_crossover, else hessenberg_unblocked. | |
| double | small_reflector (double(&v)[3], std::size_t count) |
The Householder reflection I − τ v vᵀ, v = (1, x…), that takes (α, x…) to (β, 0…), as LAPACK's dlarfg for two or three entries: on return v[0] = β and v[1…] holds x's part of v. | |
| std::complex< double > | complex_reflector (std::complex< double >(&v)[2]) |
The Householder reflection I − τ v vᴴ, v = (1, x), whose adjoint takes (α, x) to (β, 0), β real, as LAPACK's zlarfg for two entries: on return v[0] = β and v[1] holds x's part of v. | |
| standard_block | standardize_block (double &a, double &b, double &c, double &d) |
The Schur factorisation of the real 2 × 2 block [a b; c d], as LAPACK's dlanv2: [a b; c d] ← [cs −sn; sn cs]ᵀ [a b; c d] [cs −sn; sn cs], upper triangular where its eigenvalues are real, else with equal diagonal entries and b c < 0. | |
| template<typename T > | |
| void | rotate (T *x, T *y, std::size_t count, std::size_t stride, double cs, double sn) |
(x, y) ← (cs x + sn y, cs y − sn x) over count pairs stride apart: a plane rotation. | |
| template<typename T > | |
| void | reflect_three_rows (T *r1, T *r2, T *r3, std::size_t from, std::size_t to, const T &tau, const T &v2, const T &v3, bool adjoint) |
A reflection I − τ v vᴴ, v = (1, v₂, v₃), on the entries from … to − 1 of three rows: its adjoint from the left with adjoint (rows of H), else from the right (rows of Zᵀ, which are columns of Z). A reflection of two passes r3 = r2 and v3 = 0. One loop over the entries, which vectorises. | |
| template<typename T > | |
| void | reflect_two_rows (T *r1, T *r2, std::size_t from, std::size_t to, const T &tau, const T &v2, bool adjoint) |
reflect_three_rows for a reflection of two entries, v = (1, v₂), on two rows: its own loop, which vectorises where the third row aliasing the second would not. | |
| template<typename T > | |
| void | reflect_three_columns (T *first, std::size_t rows, std::size_t ld, const T &tau, const T &v2, const T &v3) |
A reflection I − τ v vᴴ, v = (1, v₂, v₃), from the right on rows rows (ld apart) of the three columns at first: x ← x − τ (x₁ + v₂ x₂ + v₃ x₃) (1, v̄₂, v̄₃) on each. The rows are short and apart: one call for them all. | |
| template<typename T > | |
| void | reflect_two_columns (T *first, std::size_t rows, std::size_t ld, const T &tau, const T &v2) |
The same for a reflection of two entries, v = (1, v₂), on two columns. | |
| void | real_schur (const schur_target< double > &t, buffer< double > &re, buffer< double > &im, std::size_t steps) |
The real Schur form of the upper Hessenberg t's active block, as LAPACK's dlahqr: Francis's implicit double-shift QR, with the deflation criterion of Ahues and Tisseur and exceptional shifts after 10 and 20 iterations without deflation. The eigenvalues are left in re and im, a complex pair's in two consecutive places, the positive imaginary part first. | |
| void | complex_schur (const schur_target< std::complex< double > > &t, buffer< std::complex< double > > &values, std::size_t steps) |
The complex Schur form of the upper Hessenberg t's active block, as LAPACK's zlahqr: the shifted QR algorithm with Wilkinson's shift, the deflation criterion of Ahues and Tisseur, and exceptional shifts after 10 and 20 iterations without deflation. The eigenvalues are left in values. | |
| std::size_t | count_of (lapack_index count) |
| A count as an unsigned size (never negative where it is taken). | |
| constexpr lapack_index | multishift_from (double) |
| The order from which a real Schur form is taken by multishift_schur, below which real_schur is faster (eight arm64 cores: order 128 1.85 against 2.60 ms; 224, 8.9 against 8.6; with the vectors, 192, 8.1 against 7.1). | |
| constexpr lapack_index | multishift_from (const std::complex< double > &) |
| The same for a complex one (order 76 1.8 against 2.1 ms, 128 7.5 against 6.5; with the vectors, 100, 6.4 against 5.6). | |
| plane_rotation< double > | givens (double f, double g) |
The rotation that takes (f, g) to (r, 0), as LAPACK's dlartg: r of f's sign. | |
| plane_rotation< std::complex< double > > | givens (const std::complex< double > &f, const std::complex< double > &g) |
The rotation that takes (f, g) to (r, 0), as LAPACK's zlartg: c real, r of f's phase. | |
| template<typename T > | |
| void | rotate_pair (T *x, T *y, std::size_t count, std::size_t stride, double c, const T &s) |
(x, y) ← (c x + s y, c y − s̄ x) over count pairs stride apart. | |
| template<typename T > | |
| void | exchange_diagonal (one_based< T > t, one_based< T > qt, lapack_index n, lapack_index k) |
Exchanges the diagonal entries k and k + 1 of the upper triangular (or quasi-triangular, both 1 × 1 there) t by a plane rotation, as LAPACK's ztrexc step and dlaexc for two 1 × 1 blocks; the rows of qt (Qᵀ) turn with it. | |
| void | reflect_rows_3 (const std::array< double, 3 > &u, double tau, double *t, std::size_t ld, std::size_t count) |
H = I − τ u uᵀ of three entries applied from the left to the 3 × count rows at t (row stride ld). | |
| void | reflect_columns_3 (const std::array< double, 3 > &u, double tau, double *t, std::size_t ld, std::size_t count) |
The same from the right, to the count × 3 columns at t. | |
| void | small_sylvester (const one_based< double > &tl, const one_based< double > &tr, const one_based< double > &b, lapack_index n1, lapack_index n2, double(&x)[4]) |
The solution X (n1 × n2, row-major in x) of TL X − X TR = B, n1 and n2 at most 2, by Gaussian elimination with complete pivoting on its Kronecker form, a pivot below ε times the largest coefficient taken as that, as LAPACK's dlasy2 (the blocks' entries near 1: no scaling against overflow). | |
| void | standardize_at (one_based< double > t, one_based< double > qt, lapack_index n, lapack_index j) |
The standard form of the 2 × 2 block at j of the quasi-triangular t, its rotation carried to the rest of T and to the rows of qt, as dlaexc ends. | |
| bool | exchange_blocks (one_based< double > t, one_based< double > qt, lapack_index n, lapack_index j1, lapack_index n1, lapack_index n2) |
Exchanges the adjacent diagonal blocks of the real quasi-triangular t at j1, of orders n1 then n2, by an orthogonal similarity carried to the rows of qt, as LAPACK's dlaexc. | |
| lapack_index | block_order (one_based< double > t, lapack_index n, lapack_index j) |
The order of the diagonal block of the quasi-triangular t that starts at j: 2 where a subdiagonal entry below it is nonzero. | |
| bool | reorder (one_based< double > t, one_based< double > qt, lapack_index n, lapack_index ifst, lapack_index &ilst) |
Moves the diagonal block of the real quasi-triangular t at ifst to ilst by exchanges of adjacent blocks, the rows of qt (Qᵀ) turning with it, as LAPACK's dtrexc. ilst becomes where the block ended: short of its goal where an exchange was refused (then. | |
| bool | reorder (one_based< std::complex< double > > t, one_based< std::complex< double > > qt, lapack_index n, lapack_index ifst, lapack_index &ilst) |
Moves the diagonal entry of the complex upper triangular t at ifst to ilst, as LAPACK's ztrexc. | |
| void | bulge_column (one_based< double > h, lapack_index count, const std::complex< double > &s1, const std::complex< double > &s2, double *v) |
A multiple of the first column of (H − s₁ I)(H − s₂ I) for the leading 2 × 2 or 3 × 3 (count) of the real Hessenberg h, the shifts a conjugate pair or two reals, as LAPACK's dlaqr1. | |
| void | bulge_column (one_based< std::complex< double > > h, lapack_index count, const std::complex< double > &s1, const std::complex< double > &s2, std::complex< double > *v) |
The same for a complex h and any two shifts, as LAPACK's zlaqr1. | |
| void | pair_shifts (std::complex< double > *s, lapack_index count, double) |
Arranges a real matrix's shifts so that each pair is a conjugate pair or two reals, as dlaqr5 begins (conjugate shifts already adjacent). | |
| void | pair_shifts (std::complex< double > *, lapack_index, const std::complex< double > &) |
| A complex matrix's shifts pair as they come. | |
| double | as_scalar (const std::complex< double > &s, double) |
The shift s as the matrix's scalar: a real matrix's are real where they are used one by one. | |
| std::complex< double > | as_scalar (const std::complex< double > &s, const std::complex< double > &) |
| A complex matrix's shift as it is. | |
| template<typename T > | |
| void | vigilant_deflation (one_based< T > h, lapack_index k, lapack_index ktop, lapack_index kbot, double smallest) |
H(k + 1, k) set to zero where the deflation criterion of Ahues and Tisseur holds, as the sweep of dlaqr5 tests after each step of a bulge (vigilant deflation). | |
| template<typename T > | |
| void | reflect_from_left (T *x, std::size_t stride, const T &tau, const T &v2, const T &v3) |
The adjoint I − τ̄ v vᴴ of a reflection of three entries from the left on one column (x three entries stride apart). | |
| template<typename T > | |
| void | chase_bulges (const schur_target< T > &target, lapack_index ktop, lapack_index kbot, std::complex< double > *s, lapack_index count, bool accumulate) |
One small-bulge multishift QR sweep over the active block ktop … kbot of target's Hessenberg matrix with the count shifts at s, as LAPACK's dlaqr5 and zlaqr5: chains of 3 × 3 bulges, two shifts each, chased down the diagonal by slabs of 4 · count / 2 columns; with accumulate, each slab's reflections are accumulated in U and applied to the rest of the matrix and to Z through the matrix product kernel, else applied to them directly (LAPACK's kacc22, which accumulates from 14 shifts). A subdiagonal entry that becomes negligible on the way is set to zero (vigilant deflation). | |
| template<typename T > | |
| void | multishift_schur (const schur_target< T > &target, std::complex< double > *w, bool recursive, lapack_index steps) |
The Schur form of target's active block by the multishift QR algorithm with aggressive early deflation, as LAPACK's dlaqr0 and zlaqr0 (dlaqr4 without recursive): each iteration deflates what it can from a window at the bottom of the active block (early_deflation) and, unless that deflated enough, sweeps the block with the window's other eigenvalues as shifts (chase_bulges). The eigenvalues are left in w (one per row; a real matrix's complex pairs adjacent, the positive imaginary part first). | |
| void | plain_schur (const schur_target< double > &t, std::complex< double > *w) |
The Schur form of t's active block by real_schur, its eigenvalues to w (one per row of t). | |
| void | plain_schur (const schur_target< std::complex< double > > &t, std::complex< double > *w) |
| The same by complex_schur. | |
| lapack_index | undeflated (one_based< double > t, one_based< double > vt, lapack_index jw, double s, double smallest, std::complex< double > *w) |
The eigenvalues of the real Schur form t (jw × jw) of a deflation window that do not deflate, moved to its top, as dlaqr3: from the bottom, a block whose spike entries s v₁ⱼ are negligible deflates, any other is moved up; the undeflated blocks are then sorted by decreasing modulus (which helps graded matrices), and every eigenvalue of the window written to w. | |
| lapack_index | undeflated (one_based< std::complex< double > > t, one_based< std::complex< double > > vt, lapack_index jw, const std::complex< double > &s, double smallest, std::complex< double > *w) |
The same for a complex Schur form, as zlaqr3: 1 × 1 blocks only, sorted by selection. | |
| template<typename T > | |
| deflation_count | early_deflation (const schur_target< T > &target, lapack_index ktop, lapack_index kbot, lapack_index nw, std::complex< double > *w, bool recursive) |
Aggressive early deflation, as LAPACK's dlaqr3 and zlaqr3: the trailing nw × nw window of the active block ktop … kbot is brought to Schur form; its eigenvalues whose coupling to the rest (the spike, s times the first row of the window's Schur vectors) is negligible deflate; the window is reduced back to Hessenberg form above them, and its transformation carried to the rest of H and to Z through the matrix product kernel. Every eigenvalue of the window is written to w: the deflated ones final, the others the next shifts. recursive: the window itself by multishift_schur past its crossover. | |
| lapack_index | recommended_shifts (lapack_index nh) |
The shifts LAPACK's iparmq recommends for an active block of order nh. | |
| lapack_index | recommended_window (lapack_index nh, lapack_index ns) |
The deflation window LAPACK's iparmq recommends for an active block of order nh and ns shifts. | |
| void | exceptional_shifts (one_based< double > h, std::complex< double > *w, lapack_index ks, lapack_index kbot, lapack_index ktop) |
Exceptional shifts for a real matrix, as dlaqr0's after kexsh iterations without deflation: from the subdiagonal entries, as conjugate pairs of the 2 × 2 blocks dlahqr's exceptional shifts build. | |
| void | exceptional_shifts (one_based< std::complex< double > > h, std::complex< double > *w, lapack_index ks, lapack_index kbot, lapack_index) |
Exceptional shifts for a complex matrix, as zlaqr0's. | |
| void | pair_for_sweep (std::complex< double > *w, lapack_index ks, lapack_index kbot, double) |
A real matrix's shifts brought into pairs, each a conjugate pair or two reals, as dlaqr0 does before its sweep (conjugate shifts already adjacent). | |
| void | pair_for_sweep (std::complex< double > *, lapack_index, lapack_index, const std::complex< double > &) |
| A complex matrix's shifts need no pairing. | |
| bool | two_real_shifts (const std::complex< double > *w, lapack_index ks, lapack_index kbot, double) |
| Whether a real matrix's sweep is down to two shifts, both real: then the one nearer the corner serves twice. | |
| bool | two_real_shifts (const std::complex< double > *, lapack_index ks, lapack_index kbot, const std::complex< double > &) |
| The same for a complex matrix: two shifts of any kind. | |
| void | schur_form (const schur_target< double > &t, buffer< double > &re, buffer< double > &im, std::size_t steps) |
The real Schur form of t's active block: real_schur below multishift_from, else multishift_schur, as LAPACK's dhseqr chooses. steps bounds real_schur's iterations. | |
| void | schur_form (const schur_target< std::complex< double > > &t, buffer< std::complex< double > > &values, std::size_t steps) |
| The complex Schur form likewise: complex_schur or multishift_schur. | |
| template<typename S > | |
| double | solve_block (const S(&m)[4], S(&r)[2], double floor) |
Solves [a b; c d] x = r in place of r by Gaussian elimination with complete pivoting, a pivot below floor taken as floor, as LAPACK's dlaln2; r is scaled first, and. | |
| template<typename S , typename P > | |
| scaled_quotient< S > | divided (S r, P pivot, double floor) |
r / pivot, a pivot below floor taken as floor, r scaled first where the quotient would pass back_substitution_limit. | |
| void | normalize_vector (std::complex< double > *v, std::size_t n, bool signed_real) |
Scales v to unit Euclidean norm and turns it so that its first component of largest modulus is real: positive for a complex matrix (zgeev), of its real part's sign for a real one's pair (dgeev). | |
| template<typename T , typename X > | |
| void | solve_upward (const buffer< T > &t, std::size_t n, std::size_t k, std::size_t top, typename X::scalar lambda, double floor, X &x) |
Solves (T − λ I) x = 0 upward from the component k, rows k − 1 … 0 of the row-major quasi-triangular t (its 2 × 2 blocks where a subdiagonal entry is nonzero), x [k … top] given: as LAPACK's dtrevc and ztrevc, a pivot below floor taken as floor, the components scaled down where they would grow past back_substitution_limit. The sums of rows against x are X 's, vectorised. | |
| buffer< std::complex< double > > | complex_vectors (const buffer< std::complex< double > > &t, std::size_t n) |
The eigenvectors of the complex upper triangular t (row-major n × n) as the rows of Xᵀ: row k is x with (T − t_kk I) x = 0, x_k = 1, zero past k, as LAPACK's ztrevc. | |
| real_eigenvectors | real_vectors (const buffer< double > &t, const buffer< double > &re, const buffer< double > &im, std::size_t n) |
The eigenvectors of the real quasi-triangular t (row-major n × n, its 2 × 2 blocks standard) as the rows of Xᵀ, as LAPACK's dtrevc: real for a real eigenvalue, solved in real arithmetic; for a pair re ± i im (im [k] > 0 the first) the first's vector, the second's its conjugate. | |
| template<typename T > | |
| buffer< T > | back_transformed (const buffer< T > &xt, const buffer< T > &zt, std::size_t n) |
Vᵀ = Xᵀ Zᵀ, row-major n × n, through the matrix product kernel. | |
| bool | gemm_path_usable (gemm_path path) noexcept |
Whether path can run here. | |
| gemm_path | best_gemm_path (gemm_path most) noexcept |
The most capable path that can run here, most at the highest. | |
| std::atomic< gemm_path > & | gemm_path_cap () noexcept |
| The most capable path products may take, as numerical::limit_gemm_path set it: none is probed until a product (or a question) reaches it, so a cap below the AMX keeps its instructions from ever running. | |
| std::atomic< std::uint8_t > & | gemm_path_found () noexcept |
The path gemm_path_current found last, packed with the cap it was found under (cap · 8 + path); 0xFF before any. | |
| gemm_path | gemm_path_current () noexcept |
The path the next double product takes: the most capable one under the cap that runs here (each path's probe is made once). Kept with its cap, one load, where each call walked the paths down from the cap through their probes' guards: every small step of a decomposition asks (an LU of order 100 through dgetrf_ on one AVX2 core spent 5 % of its time here). A race stores the same value twice. | |
| void | subtract_product (double *y, const double &f, const double *x, std::size_t count) |
y[i] −= f · x[i] for double on the widest x86-64 units the matrix product's kernel may use (gemm_path_current), each element one fused multiply-add there; elsewhere subtract_product's template. The decompositions' row updates take it: a wheel built for the baseline x86-64 holds two double per vector and fuses none. | |
| void | subtract_product (double *y, const double *x, const double &f, std::size_t count) |
y[i] −= x[i] · f for double, as the other (a real product is the same either way round). | |
| void | subtract_products (double *y, const double *f, const double *x, std::size_t ldx, std::size_t terms, std::size_t count) |
y[i] -= Σₜ f[t] · x[(t · ldx) + i] for double, the units chosen as subtract_product chooses them, the elements held in registers through the terms there. | |
| void | subtract_rows (double *y, std::size_t ldy, const double *f, std::size_t ldf, const double *x, std::size_t ldx, std::size_t terms, std::size_t count) |
subtract_rows for double, the units chosen as subtract_products chooses them: each element the value subtract_products gives it. | |
| void | subtract_products (std::complex< double > *y, const std::complex< double > *f, const std::complex< double > *x, std::size_t ldx, std::size_t terms, std::size_t count) |
subtract_products for std::complex<double>: on AVX-512 or AVX2 with FMA where the matrix product's kernel runs there, each part of each term one fused multiply-add (Re f · x, then Im f times the swapped parts), the elements in registers through the terms; elsewhere the template, whole complex products as the compiler rounds them. | |
| void | subtract_product (std::complex< double > *y, const std::complex< double > &f, const std::complex< double > *x, std::size_t count) |
y[i] −= f · x[i] for std::complex<double>: one term of subtract_products on AVX-512 or AVX2 where the matrix product's kernel runs there, else the template. | |
| void | subtract_product (std::complex< double > *y, const std::complex< double > *x, const std::complex< double > &f, std::size_t count) |
y[i] −= x[i] · f for std::complex<double>, as the other (the template keeps the expression's order). | |
| void | add_product (double *y, const double &f, const double *x, std::size_t count) |
y[i] += f · x[i] for double, as subtract_product. | |
| void | add_product (std::complex< float > *y, const std::complex< float > &f, const std::complex< float > *x, std::size_t count) |
y[i] += f · x[i] for complex float: on the x86-64 paths by AVX2 vectors of its arithmetic (avx2_scaled_add_c32), the same bits as add_product's template, which runs elsewhere. | |
| void | add_product (std::complex< double > *y, const std::complex< double > &f, const std::complex< double > *x, std::size_t count) |
add_product for complex double (avx2_scaled_add_c64 on the x86-64 paths). | |
| double | dot_product (const double *x, const double *y, std::size_t count) |
Σ x[i] · y[i] for double, as subtract_product chooses its units (in another order of summation there, as fixed). | |
| template<split_complex C> | |
| C | complex_dot_product (const C *x, const C *y, std::size_t count) |
Σ x[i] · y[i] of complex numbers: on the x86-64 units with AVX2 (the AVX-512 path too) their fused chains (avx2_complex_dot), elsewhere the native vectors' (in another order of summation there, as fixed). | |
| std::complex< double > | dot_product (const std::complex< double > *x, const std::complex< double > *y, std::size_t count) |
dot_product for std::complex<double>, through complex_dot_product. | |
| std::complex< float > | dot_product (const std::complex< float > *x, const std::complex< float > *y, std::size_t count) |
dot_product for std::complex<float>, through complex_dot_product. | |
| template<split_complex C> | |
| void | complex_row_dots (const C *m, std::size_t ld, std::size_t rows, const C *x, std::size_t count, C *out) |
| row_dots of complex numbers: each row by complex_dot_product, as chunked_dot sums a row. | |
| void | row_dots (const std::complex< double > *m, std::size_t ld, std::size_t rows, const std::complex< double > *x, std::size_t count, std::complex< double > *out) |
row_dots for std::complex<double>, through complex_row_dots. | |
| void | row_dots (const std::complex< float > *m, std::size_t ld, std::size_t rows, const std::complex< float > *x, std::size_t count, std::complex< float > *out) |
row_dots for std::complex<float>, through complex_row_dots. | |
| void | row_dots (const double *m, std::size_t ld, std::size_t rows, const double *x, std::size_t count, double *out) |
row_dots for double: on the x86-64 units each row by their dot_product (wider than four rows of the portable vectors), elsewhere four rows at a time. | |
| template<typename T > | |
| void | gemm_kernel (std::size_t kc, const T *ap, const T *bp, std::size_t ldb, gemm_tile< T > &tile) |
tile += Aₚ·Bₚ over one panel, for any scalar. | |
| template<std::size_t Mr, typename T > | |
| void | gemm_pack_a (const T *a, std::size_t lda, std::size_t rows, std::size_t kc, T *out) |
Copies rows × kc of A (from row-major, leading dimension lda) into slivers of Mr rows, each stored column by column; a last partial sliver is padded with zeros. | |
| template<std::size_t Nr, typename T > | |
| void | gemm_pack_b (const T *b, std::size_t ldb, std::size_t kc, std::size_t cols, T *out) |
Copies kc × cols of B (from row-major, leading dimension ldb) into slivers of Nr columns, each stored row by row; a last partial sliver is padded with zeros. | |
| template<std::size_t Nr, typename T > | |
| void | gemm_pack_bt (const T *bt, std::size_t ldbt, std::size_t kc, std::size_t cols, T *out) |
gemm_pack_b of B = Bᵀᵀ, from the row-major bt (cols × kc, leading dimension ldbt): each column of a sliver is a row of bt, read along it; a last partial sliver is padded with zeros. | |
| template<typename T > | |
| const T * | rows_in_place (const T *x, std::size_t ld, std::size_t &ld_out) noexcept |
A row-major rows × cols operand for the plain loop, x itself (ld apart). | |
| template<typename T > | |
| const T * | rows_from_transpose (const T *x, std::size_t ld, std::size_t rows, std::size_t cols, buffer< T > &scratch, std::size_t &ld_out) |
A row-major rows × cols operand for the plain loop, copied into scratch from its transpose x (cols × rows, ld apart). | |
| template<std::size_t Mr, typename T > | |
| void | gemm_pack_at (const T *at, std::size_t ldat, std::size_t rows, std::size_t kc, T *out) |
The slivers of Mr rows of A (rows × kc), column by column, from its transpose at (kc × rows row-major, ldat apart): each column of a sliver is a run of a row of at; a last partial sliver is padded with zeros, as gemm_pack_a lays them out. | |
| template<typename T > | |
| T | rounded_apart_add (T sum, T a, T b) noexcept |
sum + a·b as the vector tile rounds it (vector_tile_kernel), so that the plain loop's elements are its own: under Clang's -ffp-contract=on one expression, fused in both; under GCC's, whose vector extensions fuse nothing while a scalar statement's multiply-add fuses, the product apart, a statement of its own. | |
| template<typename R > | |
| std::complex< R > | rounded_apart_add (std::complex< R > sum, std::complex< R > a, std::complex< R > b) noexcept |
| rounded_apart_add of complex numbers: the product in line (operator_times, the operator's own rounding, unfused), where it was a call that took a sixth of a complex LU. | |
| template<typename T > | |
| T | rounded_term (rounding_is< term_rounding::compiler >, T sum, T a, T b) noexcept |
sum + a·b as the compiler rounds it (rounded_apart_add). | |
| template<typename T > | |
| T | rounded_term (rounding_is< term_rounding::fused >, T sum, T a, T b) noexcept |
sum + a·b as one fused multiply-add. | |
| template<std::size_t R, std::size_t N, term_rounding Mode, typename T > | |
| void | narrow_block (std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) noexcept |
C += A·B for R rows and N columns of C, each sum in ascending order of the depth, as the plain loop's: the block held in registers, of a size the compiler knows; each term rounded as Mode says. Out of line: its unrolled sizes, inlined, spent GCC's growth budget, and the LU factorisation, which runs none of them, lost 3 to 8 % at every order (one arm64 thread, GCC 14). | |
| template<std::size_t R, term_rounding Mode, typename T > | |
| void | narrow_rows (std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) noexcept |
C += A·B for R whole rows of C of n columns: blocks of narrow_columns, the last narrower. | |
| template<term_rounding Mode, typename T > | |
| void | gemm_rows_rounded (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
C += A·B by the plain i-k-j loop, parallel over the rows of C. Past a depth of 1 (narrow_least_depth), the rows run narrow_group at a time by blocks (narrow_rows), the last rows one by one: the same sums, six to ten times as fast from order 4 to 16 (M1 Pro) as the loop of unknown width, whose sums the compiler cannot keep in registers. A depth of 1 streams as fast through that loop. | |
| bool | rows_fuse () noexcept |
Whether the tile a double or float product past the plain loop runs on fuses each term (gemm_rows rounds its own so): double on arm64's NEON and AMX paths, float on the NEON tile there. The x86-64 tiles' double takes term_rounding::x86_fused; the vector tile, and so the portable path, rounds as the compiler does, as the plain loop then does. On arm64 every cap but the portable one runs the NEON or the AMX path, so the cap alone answers: one load, where the path in use took an order-4 product on one arm64 core from 10.6 to 13.3 ns (the cap, 11.5). | |
| template<typename T > | |
| void | gemm_rows (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
C += A·B by the plain i-k-j loop (gemm_rows_rounded), each element rounded as the tile of the path in use rounds it, so that it is the same whatever the shape of the product it falls in: a caller that splits its columns by the threads (an LU inverse) then has the same bits on any number of them. | |
| template<typename T > | |
| T * | line_aligned (T *at) noexcept |
The first element of at's run that starts a cache line (or at's nearest element towards one, for an element whose size does not divide the allocator's alignment). | |
| template<typename Kernel > | |
| constexpr bool | gemm_shares (std::size_t m, std::size_t n, std::size_t k) noexcept |
Whether the product's threads pay: below Kernel's serial_below multiply-adds, the product runs on the caller alone. | |
| template<typename Kernel , typename T > | |
| void | prefetch_next_row (const T *, std::size_t) noexcept |
| Nothing, for a tile that prefetches no C (prefetches_next_row_of_c). | |
| template<typename Kernel , typename T > requires prefetches_next_row_of_c<Kernel> | |
| void | prefetch_next_row (const T *c, std::size_t ldc) noexcept |
The C tile at c prefetched by Kernel's prefetch_rows. | |
| template<typename Kernel > | |
| constexpr std::size_t | task_rows_cap (std::size_t) noexcept |
The most rows of a task of Kernel over a panel of depth kc: its block of rows. | |
| template<typename Kernel > requires requires(std::size_t kc) { Kernel::task_rows_for(kc); } | |
| constexpr std::size_t | task_rows_cap (std::size_t kc) noexcept |
task_rows_cap for a tile whose task rows depend on the panel's depth (Kernel::task_rows_for). | |
| template<typename Kernel , typename T , typename B = b_rows, typename A = a_rows> | |
| void | gemm_packed (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc, const T *a_packed=nullptr) |
The packed product with the tile Kernel computes: for each block of B's columns and each panel of the depth, both operands are packed (in parallel, by the blocks the product takes), then tasks of gemm_mc rows by gemm_slivers_per_task slivers of columns run in parallel, each writing a disjoint block of C, one tile at a time. Given a_packed (A as pack_left packs it for Kernel, a depth of one panel), A is read from there and only B is packed: the same tiles, the same bits. B is B's layout (b_rows, b_transposed), A A's (a_rows, a_transposed): the packing alone differs, so the bits do not. | |
| template<typename Kernel , typename T , typename B = b_rows, typename A = a_rows> | |
| void | gemm_packed_triangle (std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc, bool lower) |
gemm_packed restricted to a triangle of the square C = A·B (n × n): its lower one (j ≤ i) where lower, else its upper one (j ≥ i); the other is neither written nor depended on. The same packing (A's row blocks with no tile in the triangle left out) and the same tiles, so each element's sum is the product's; a tile across the diagonal runs on a copy of its block of C, the triangle's part copied back. A syrk takes the triangle in one pass, its operands packed once a panel, where products by halves packed them again for every block. | |
| template<typename T > | |
| bool | gemm_elsewhere (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
| A product a faster path than the vector one takes, if any: none for this scalar. | |
| template<typename B , typename A = a_rows> | |
| bool | gemm_double_elsewhere (std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, const double *b, std::size_t ldb, double *c, std::size_t ldc) |
A double product on the path in use (numerical::gemm_path_in_use), unless it is the portable one, B and A laid out as B and A say (b_rows, a_rows and their transposes). | |
| template<typename A , typename B > | |
| void | gemm_triangle (std::size_t n, std::size_t k, const double *a, std::size_t lda, const double *b, std::size_t ldb, double *c, std::size_t ldc, bool lower) |
The triangle of the double product C = A·B (gemm_packed_triangle) on the tile of the path in use, A and B laid out as A and B say. | |
| bool | gemm_elsewhere (std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, const double *b, std::size_t ldb, double *c, std::size_t ldc) |
| gemm_double_elsewhere of a row-major B. | |
| template<typename B , typename A = a_rows> | |
| bool | gemm_float_elsewhere (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc) |
A float product on the AMX where it suits the product, on the NEON tile on arm64's fused paths, on the x86-64 tiles there; elsewhere the vector tile. B and A laid out as B and A say. | |
| bool | gemm_elsewhere (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc) |
| gemm_float_elsewhere of row-major operands. | |
| template<typename T > | |
| void | gemm (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
C += A·B, every element of C summed in ascending order of the depth. | |
| bool | split_by_parts (std::size_t m) noexcept |
Whether gemm_split runs a product of m rows of A by four real products of each block's parts (split_parts_block), else by one of its expanded operands (split_expanded_block): the same bits, the parts' products faster where they read B's parts in place (gemm_packed). On the AMX path up to split_parts_rows; elsewhere (the NEON path) up to one row of the NEON tile, its b_in_place_rows: M1 Pro, all threads, by 400 × 500, 1 to 4 rows 2 to 18 % faster, 6 to 32 rows 5 to 52 % slower. | |
| template<typename R > | |
| split_operands< R > & | split_storage () |
| The calling thread's split_operands. A product uses them whole on one thread, and the real product it runs on them nests no complex one. | |
| template<split_complex Z, typename B = b_rows, typename A = a_rows> | |
| void | split_parts_block (std::size_t m, std::size_t nb, std::size_t kb, std::size_t pc, std::size_t jc, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, typename Z::value_type *ar, typename Z::value_type *br, typename Z::value_type *cr) |
One block of gemm_split ([Re C | Im C] at cr, m × nb each, += A·B, A m × kb and B kb × nb) by four real products of its parts: Re C += Re A·Re B, then −Im A·Im B; Im C += Re A·Im B, then Im A·Re B, each continuing its elements' chains ((−x)·(−y) is x·y exactly: the bits of −Im A·−Re B). B's parts are a copy of its size, which a product of a few rows of A reads in place. The block is A's columns and B's rows from pc, B's columns from jc, A and B laid out as A and B say. | |
| template<split_complex Z, typename B = b_rows, typename A = a_rows> | |
| void | split_expanded_block (std::size_t m, std::size_t nb, std::size_t kb, std::size_t pc, std::size_t jc, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, typename Z::value_type *ar, typename Z::value_type *br, typename Z::value_type *cr) |
One block of gemm_split (split_parts_block's) by one real product of its expanded operands, [Re C | Im C] += [Re A | −Im A] · [[Re B, Im B], [Im B, −Re B]]: the same chains, B copied four times its size, the product as wide and as deep as the block's two. | |
| template<split_complex Z, typename B = b_rows, typename A = a_rows> | |
| void | gemm_split (std::size_t m, std::size_t n, std::size_t k, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, Z *c, std::size_t ldc) |
C += A·B for complex numbers through the real kernel, by blocks of split_columns columns and split_depth of depth, on the real path's units and threads (a scalar complex tile has none). Within a block, each element of C takes the products of the real parts of A in ascending order of the depth, then those of the imaginary parts: by four real products of the block's parts (split_parts_block) or one of its expanded operands (split_expanded_block), as split_by_parts chooses, the same bits either way. A and B laid out as A and B say, read where they lie as the parts are copied. | |
| bool | plain_loop_pays (std::size_t m, std::size_t n, std::size_t k) noexcept |
Whether a double product of m rows, n columns and depth k, past gemm_rows's threshold, still runs faster in the plain loop than packed on the path in use: packing B's slivers for so few rows of A, or so short a depth, did not pay. One core, the plain loop's time over the packed product's: the NEON tile at a depth of 16 or less whatever the shape (0.5 to 0.6 up to 2000 × 2000), at 64 or less for 64 rows or fewer (0.4 to 0.9; the same under GCC 14 and on eight threads); the AVX2 tile for 16 rows or fewer to a depth of 64, or 32 rows to 32 (0.2 to 0.8, six cores too), as wide as plain_loop_least_width. The AMX's packing pays from 16 rows. | |
| template<typename T > | |
| bool | plain_loop_suits (std::size_t, std::size_t, std::size_t) noexcept |
plain_loop_pays for a product of T: double's alone, the one it was measured for. | |
| template<> | |
| bool | plain_loop_suits< double > (std::size_t m, std::size_t n, std::size_t k) noexcept |
plain_loop_pays for double. | |
| template<typename T > | |
| void | gemm_large (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
| A product past the plain loop: on the AMX or x86-64 tiles where gemm_elsewhere takes it, else on the vector tile. | |
| template<typename B , typename A , split_complex T> | |
| void | gemm_large_complex (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
A complex product past the plain loop: on the vector tile on the portable path, whole complex products summed in order; on the x86-64 complex tiles there; else (arm64) through the real kernel (gemm_split). B and A laid out as B and A say. | |
| template<typename B , typename A > | |
| void | gemm_large_layouts (std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, const double *b, std::size_t ldb, double *c, std::size_t ldc) |
A double product past the plain loop, B and A laid out as B and A say (gemm_layouts). | |
| template<typename B , typename A > | |
| void | gemm_large_layouts (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc) |
A float product past the plain loop, B and A laid out as B and A say (gemm_layouts). | |
| template<typename B , typename A , split_complex T> | |
| void | gemm_large_layouts (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
A complex product past the plain loop, B and A laid out as B and A say (gemm_layouts). | |
| template<typename A , typename B , typename T > | |
| void | gemm_layouts (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
gemm with either operand given as its transpose (A: a_rows or a_transposed, B: b_rows or b_transposed), so that a product by a transpose (a syrk's P Pᵀ, BLAS's transa and transb, a NumPy array in Fortran order) needs no copy of it. The tiles are gemm's, the operands packed (or, through gemm_split, copied into their parts) from where they lie, and so are the bits; a product gemm runs in the plain loop copies a transposed operand row-major and runs there, small as it is. | |
| template<typename T > | |
| void | gemm_bt (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *bt, std::size_t ldbt, T *c, std::size_t ldc) |
gemm_layouts with B given as Bᵀ (bt, n × k, ldbt apart). | |
| template<typename T > | |
| void | gemm_at (std::size_t m, std::size_t n, std::size_t k, const T *at, std::size_t ldat, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
gemm_layouts with A given as Aᵀ (at, k × m, ldat apart). | |
| template<typename T > | |
| void | gemm_atbt (std::size_t m, std::size_t n, std::size_t k, const T *at, std::size_t ldat, const T *bt, std::size_t ldbt, T *c, std::size_t ldc) |
| gemm_layouts with both operands given as their transposes. | |
| template<typename T > | |
| packed_left< T > | read_in_place (std::size_t m, std::size_t k, const T *a, std::size_t lda) noexcept |
The m × k a read in place by gemm_with, unpacked: each product packs it as gemm does. | |
| template<typename T > | |
| packed_left< T > | left_rows (const packed_left< T > &left, std::size_t r0, std::size_t rows) noexcept |
The rows [r0, r0 + rows) of left, r0 a multiple of left.granule: their slivers in its packing, or their rows of A in place. | |
| template<typename T > | |
| left_kernel | left_kernel_for (std::size_t, std::size_t, std::size_t) noexcept |
The tile gemm runs a product of m × n by depth k on, for a left operand packed once for many products: plain for a scalar other than double (written as an overload per scalar, as conj_of). | |
| template<> | |
| left_kernel | left_kernel_for< double > (std::size_t m, std::size_t n, std::size_t k) noexcept |
left_kernel_for a double product: past gemm_rows, as gemm_elsewhere chooses, of a depth of one panel; plain otherwise. | |
| template<typename Kernel , typename T > | |
| const T * | pack_left_for (std::size_t m, std::size_t k, const T *a, std::size_t lda, buffer< T > &storage) |
A packed for Kernel into storage, in parallel by its blocks of rows, as gemm_packed packs it. | |
| template<typename T > | |
| packed_left< T > | pack_left (std::size_t m, std::size_t, std::size_t k, const T *a, std::size_t lda, buffer< T > &) |
The m × k a packed once for products with it of about n columns each (gemm_with), in storage (grown as needed, kept by the caller from one call to the next): on the tile gemm would choose for m × n × k. A scalar other than double (here), a depth of more than one panel, or a product gemm_rows takes, leaves A to be read in place. | |
| template<> | |
| packed_left< double > | pack_left< double > (std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, buffer< double > &storage) |
pack_left for double, on the tile left_kernel_for chooses. | |
| template<typename T > | |
| void | gemm_with (const packed_left< T > &left, std::size_t n, const T *b, std::size_t ldb, T *c, std::size_t ldc) |
C += A·B for the packed left (pack_left) and the k × n b, as gemm computes it, A read in place here, for a scalar other than double. | |
| template<> | |
| void | gemm_with< double > (const packed_left< double > &left, std::size_t n, const double *b, std::size_t ldb, double *c, std::size_t ldc) |
gemm_with for double: on the tile A was packed for, A read from its packing, where gemm would choose that tile for this n; else as gemm. Either way gemm's tiles and order of summation, its bits. | |
| template<typename T > | |
| T | chunked_dot (const T *a, const T *b, std::size_t count) |
Σ a[i] · b[i] for a plain dot product: chunks of dot_chunk elements in parallel, each by dot_product (several vector accumulators, on the widest x86-64 units for double), their sums added in the chunks' order. The same on any number of threads; in the last bits, not the ascending sum, which waits on each addition before the next (ten times slower at 100 000 elements). | |
| template<typename T > | |
| void | matrix_times_vector (const T *a, std::size_t rows, std::size_t cols, const T *x, T *y) |
y = A x for the row-major rows × cols a: each element the dot product of its row with x, as chunked_dot sums it (a row of a chunk or less is one dot_product), the rows in parallel by groups (row_dots, which share each load of x). The same on any number of threads. A double matrix of amx_dot_least elements to fewer than amx_dot_most, on the AMX path, runs on the AMX (amx_dot_rows_f64): each element then its eight lanes of fused multiply-adds, added in order. | |
| long | current_process () noexcept |
This process's identifier: a pool started before a fork belongs to the parent, whose threads the child does not have. | |
| std::vector< std::size_t > | worker_cores (std::size_t cores, std::size_t threads) |
The core of each worker the pool starts beside its callers, from the cores known (0: none) and its threads (at least 1): one worker per core, pinned there, where the cores are known (each run then leaves out the one on the caller's core, sits_out); else threads − 1 workers free to move. | |
| constexpr bool | sits_out (std::size_t worker, std::size_t caller) noexcept |
Whether the worker pinned to core worker sits out a run whose caller is on core caller: one thread per core, as two on one core's hyperthreads share its multiply-add units. A free worker never does. | |
| std::size_t | physical_cores () noexcept |
The physical cores this process may run on, or 0 where the system does not say. Linux: core_sets; macOS: hw.physicalcpu. | |
| constexpr std::size_t | pool_threads (std::size_t logical, std::size_t physical) noexcept |
The pool's threads from the logical hardware threads and the physical cores (0: unknown): the cores where known, at least 1. | |
| std::size_t | hardware_threads () noexcept |
| The threads the pool runs, the caller's included (at least 1), known without starting the pool: one per physical core where the system tells them (else one per hardware thread). Two threads on one core's hyperthreads share its multiply-add units: a product of order 128 on eight AVX2 cores with sixteen hyperthreads took 140 µs on sixteen threads, 71 on eight; order 2048, the same either way. Asked once, as the pool's workers are started once (glibc reads a file for the hardware threads on every call). | |
| template<typename Ready > | |
| void | spin_until (const Ready &ready) |
Polls ready, yielding the core between polls, until it holds or spin_time has passed. | |
| constexpr bool | wants (std::size_t core, std::size_t caller, std::size_t helpers) noexcept |
Whether the worker on core joins a run of helpers workers whose caller is on caller: by rank among the cores, the caller's aside, so a run of few threads takes the same cores from one run to the next. | |
| std::atomic< std::size_t > & | thread_cap () noexcept |
| The most threads a parallel_for uses (0: as many as the pool has), set by numerical::set_max_threads. | |
| std::size_t | threads_for (std::size_t count, std::size_t work_per_index, std::size_t hardware) |
How many threads count indices of work_per_index each pay for: one per work_per_thread, at most one per index and per hardware thread. | |
| template<typename Body > | |
| void | parallel_for (std::size_t count, std::size_t work_per_index, const Body &body) |
Runs body(i) for every i in [0, count). | |
| template<typename Body , typename Side > | |
| void | parallel_for_beside (std::size_t count, std::size_t work_per_index, const Body &body, const Side &side) |
Runs side on the caller while the pool's workers run body(i) for every i in [0, count), then lets the caller take the chunks left: one task beside a parallel loop, each touching what the other does not (a factorisation's next panel beside its trailing update). side's own parallel loops find the pool busy and run in its thread. Where the work pays for fewer than two threads, or the pool is busy, side runs, then the loop, on the caller. An exception of side comes back after the loop's chunks have joined, before theirs. | |
| template<typename T > | |
| constexpr std::size_t | fault_chunk (std::size_t count) noexcept |
The elements one index of parallel_zero or parallel_copy takes of count: zero_chunk, or, from the size at which a buffer takes Linux's 2 MiB pages (16 MiB), one such page's worth, so that no two threads fault in the same page: they waited on each other's fault there (six x86-64 AVX2 cores, 32 MiB zeroed into fresh memory 3.6 → 2.5 ms; the LU and Cholesky of order 2000, their operand copied by eight 4 KiB pages a thread, took 2 to 4 % longer than with one thread's copy). | |
| template<typename T > | |
| void | parallel_zero (T *out, std::size_t count) |
Sets the count elements at out to zero, split across the threads the work pays for: a fresh result's pages are first touched, and faulted in, by several threads, where a sequential fill pays every fault on one (a product on eight arm64 threads: 16 % faster at order 512, 7 % at 1024). The zero is a constant in each chunk's fill, which the compiler turns into the C library's zeroing (on arm64 it clears whole cache lines without reading them; a value passed in, stored element by element, was 4 % slower at order 1024 on one thread). | |
| template<typename T > | |
| void | parallel_copy (const T *in, std::size_t count, T *out) |
Copies the count elements at in to out, split across the threads the work pays for, as parallel_zero does: a fresh copy's pages are faulted in by several threads, where a sequential copy pays every fault on one (a 2 MiB copy into fresh memory took 0.85 ms on one x86-64 core in a container, a third of the Cholesky factorisation of order 512). | |
| template<typename T > | |
| std::array< T, symmetric_row_group > | symmetric_rows_tail (const T *rows, std::size_t n, const T *v, T *out, std::size_t i, std::size_t done, std::array< T, symmetric_row_group > sums) |
The columns done…i+r of symmetric_rows, one by one: for each row r in order, out[k] += row_r[k] · v[i+r] for done ≤ k < i + r and row_r[k] · v[k] added to sums[r], its diagonal term last. | |
| template<std::size_t Rows, std::size_t Cols, typename T , typename Full > | |
| void | through_edge (T *c, std::size_t ldc, std::size_t rows, std::size_t cols, Full full) noexcept |
One tile of a matrix product kernel, rows × cols of C (at most Rows × Cols, leading dimension ldc), through full, which computes a whole tile in place given its first element and leading dimension: a whole tile in C itself, a partial one in a copy whose outside is zero. Only that copy is zeroed: zeroing it for every tile, whole ones included, took a twelfth of an x86-64 product's time. | |
| template<typename T > | |
| std::size_t | subtract_scaled_vectors (T *, const T *, const T &, std::size_t) |
| No leading elements in vectors, for a type no native vector holds. | |
| template<typename T , typename Magnitude > | |
| std::size_t | largest_magnitude_vectors (const T *, std::size_t, Magnitude &, std::size_t &) |
Only x[0], the one best and at come in with, for a type no native vector holds. | |
| template<typename T > | |
| std::size_t | subtract_combined_vectors (T *, const T *, const T *, std::size_t, std::size_t, std::size_t) |
| No leading blocks in vectors, for a type no native vector holds or a compiler that vectorises the scalar block itself (Clang). | |
| template<typename T > | |
| std::size_t | subtract_rows_vectors (T *, std::size_t, const T *, std::size_t, const T *, std::size_t, std::size_t, std::size_t) |
| No leading columns in vectors, for a type no native vector holds or a compiler that vectorises the scalar block itself (Clang). | |
| template<typename T > | |
| std::size_t | add_scaled_vectors (T *, const T *, const T &, std::size_t) |
| No leading elements in vectors, for a type no native vector holds. | |
| template<typename T > | |
| T | symmetric_row (const T *row, const T *v, const T &f, T *out, std::size_t count) |
Σ row[i] · v[i] and out[i] += row[i] · f, in order, for a type no native vector holds. | |
| template<typename T > | |
| std::array< T, symmetric_row_group > | symmetric_rows (const T *rows, std::size_t n, const T *v, T *out, std::size_t i) |
| symmetric_rows one column at a time, for a type no native vector holds. | |
| template<typename T > | |
| T | dot_product (const T *x, const T *y, std::size_t count) |
Σ x[i] · y[i], in order, for a type no native vector holds. | |
| template<typename T > | |
| T | conjugate_dot_product (const T *x, const T *y, std::size_t count) |
Σ conj(x[i]) · y[i], in order, where no native vector holds the parts. | |
| template<typename T > | |
| void | row_dots (const T *m, std::size_t ld, std::size_t rows, const T *x, std::size_t count, T *out) |
dot_product of each row of m with x, for a type no native vector holds. | |
| template<typename T > | |
| T | operator_times (const T &a, const T &b) |
a · b, as std::operator* computes it, in line. libc++'s complex product (the four products and their sums as statements of their own, so unfused) is copied here, its rescue of a product whose parts are both NaN left to it: where it was a call, out of line in a translation unit as large as the Python extension, it took 45 % of a complex LU of order 48. Elsewhere (libstdc++, whose product the compiler expands) and for a real T, the operator itself. | |
| template<typename T > | |
| void | subtract_product (T *y, const T &f, const T *x, std::size_t count) |
y[i] -= f · x[i] for count elements of y and x, which do not overlap. Each element is computed exactly as that scalar statement computes it (its product rounded before its difference, or fused with it where the compiler contracts the statement); the arguments follow the expression, since a compiler may contract f · x and x · f of complex numbers differently. Under GCC, over native vectors where T has them: at -O2 GCC leaves the scalar loop scalar, unable to rule out an overlap without a check it will not add. Clang vectorises the scalar loop itself, with that check, and better than these vectors (an LU solve of order 128 on arm64 took 19 % longer through them). For double, gemm.hpp's overload chooses the x86-64 units at run time, a fused multiply-add per element there. | |
| template<typename T > | |
| void | subtract_products (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count) |
y[i] -= Σₜ f[t] · x[(t · ldx) + i] for count elements of y, terms rows of x apart by ldx: subtract_product for each term in ascending order, so each element is computed as that sequence of calls computes it. For double, gemm.hpp's overload holds the elements in registers through the terms. | |
| template<typename T > requires std::is_floating_point_v<T> | |
| void | subtract_products (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count) |
subtract_products for a real T: eight elements of y held in a local array through the terms, which the compiler keeps in registers, where each term was a call through the whole row. Each element takes its terms in the same order as a sequence of subtract_product gives them (an LU solve of order 64 and four columns on arm64: 18 → 8.5 µs), and the same value under Clang, and under GCC on x86-64, which fuses neither; under GCC on arm64 the vectors of subtract_combined_vectors and subtract_scaled_vectors round apart and the scalar statements fused, and their splits differ. | |
| template<typename T > requires std::is_floating_point_v<T> | |
| void | subtract_rows (T *y, std::size_t ldy, const T *f, std::size_t ldf, const T *x, std::size_t ldx, std::size_t terms, std::size_t count) |
subtract_products for product_rows rows of y at once, ldy apart, each with its own factors (ldf apart in f) on the same rows of x: y_r[i] -= Σₜ f_r[t] · x[(t · ldx) + i]. Each element takes the statement of subtract_products, a term at a time in ascending order, so the same value; each element of x read serves every row, and the rows' chains of multiply-adds advance together, where a row alone waits on its last one through each term. | |
| template<typename T > | |
| void | subtract_product (T *y, const T *x, const T &f, std::size_t count) |
y[i] -= x[i] · f, as the other subtract_product (a real product is the same either way round). | |
| template<split_complex C> | |
| void | complex_scaled_add (C &y, const C &f, const C &x) |
y += f · x of one complex number, its part products rounded apart as libc++'s operator computes them (add_scaled_vectors's arithmetic, for the elements past its vectors). | |
| template<split_complex C> | |
| void | add_product (C *y, const C &f, const C *x, std::size_t count) |
| add_product of complex numbers: vectors of their parts, then complex_scaled_add, every element the same arithmetic. | |
| template<typename T > | |
| void | add_one_product (T &y, const T &f, const T &x) |
y += f · x of one element, as add_product computes each of its elements. | |
| template<split_complex C> | |
| void | add_one_product (C &y, const C &f, const C &x) |
| add_one_product of complex numbers: complex_scaled_add, as add_product computes them. | |
| template<typename T > | |
| void | add_product (T *y, const T &f, const T *x, std::size_t count) |
y[i] += f · x[i], as subtract_product. | |
| template<typename T > | |
| T | pairwise_sum (const T *first, std::size_t count, std::size_t stride) |
Sum of count values from first at stride stride, by NumPy's pairwise summation: below 8 values in order; up to 128 in eight interleaved partial sums, added pairwise, then the remainder in order; above, the two halves (the first a multiple of 8 long) summed apart. | |
| template<typename T > | |
| void | pairwise_sum_rows (const T *first, std::size_t count, std::size_t stride, std::size_t width, T *out) |
pairwise_sum of width neighbouring columns at once, into out: column q's sum that of pairwise_sum(first + q, count, stride), to the bit, each step adding a whole row of the columns (contiguous, stride apart) where a column was walked down stride at a time. | |
| template<typename T > | |
| void | pairwise_sums (const T *data, std::size_t outer, std::size_t count, std::size_t inner, T *out) |
The pairwise sums of outer blocks of count rows of inner elements (row-major), each column of a block summed down its rows, into out (outer × inner): one column at a time where inner is 1, else by rows (pairwise_sum_rows), the same bits. | |
| bool | is_label (char ch) |
Is ch an ASCII letter (a valid index label)? | |
| std::vector< std::string > | split_inputs (const std::string &lhs) |
| Splits an einsum left-hand side into per-operand subscripts. | |
| template<typename T > | |
| tensor< T > | tile_gemm (const T *a, const T *b, std::size_t rows, std::size_t inner, std::size_t cols, std::vector< std::size_t > out_shape) |
Matrix product C(M×N) = A(M×K) · B(K×N), row-major. | |
| template<typename T > | |
| tensor< T > | canonical_product (const T *a, const T *b, std::size_t rows, std::size_t inner, std::size_t cols, std::vector< std::size_t > out_shape) |
The product of the row-major rows × inner a and inner × cols b, shaped out_shape: the canonical two-operand contraction of einsum_impl (a matrix times a vector, a vector times a matrix, a matrix product, a dot product), on any storage that holds the operands contiguously. | |
| template<typename T > | |
| buffer< T > | vector_times_transposed (const T *a, const T *bt, std::size_t ldbt, std::size_t inner, std::size_t cols) |
canonical_product's vector a (inner) times a matrix given as its transpose bt (cols × inner, its rows ldbt apart, B's columns): each element of the product the sum canonical_product forms, a[k] · B[k][j] added in ascending order of k with add_product's arithmetic (add_one_product), so the same bits; a column of B a run of bt, transposed_columns of them at once. | |
| template<typename T > | |
| std::optional< tensor< T > > | sum_over_run (const std::string &in, const std::string &output, const std::vector< std::size_t > &shape, const T *data) |
einsum of one row-major operand data (shape, its axes labelled in, each label once) into the labels output, where those it sums are one run of its axes and output the others in order (ij->i, ij->, ij->j, ijk->ik): each output element the sum of evenly spaced elements, by pairwise halving as sum() adds them (pairwise_sums), where einsum's general loop added them one after another. None for any other spec (nothing summed, the summed labels apart, the output reordered). | |
| template<typename T > | |
| buffer< T > | permuted (const buffer< T > &data, const std::string &from, const std::string &to, const std::map< char, std::size_t > &sizes) |
The elements of data, whose axes carry the labels from (extents in sizes), laid out row-major along the labels to, a permutation of from that differs from it (so at least one label): an odometer over to's outer labels, a strided sweep of its last. | |
| template<typename T > | |
| std::optional< tensor< T > > | contract_pair (const std::string &sa, const tensor< T > &a, const std::string &sb, const tensor< T > &b, const std::string &output, const std::string &summed, const std::map< char, std::size_t > &sizes) |
A two-operand contraction as matrix products: the labels both operands and the output carry (batch), those one operand and the output carry (free), and those both operands carry alone (summed, in ascending order, the order the general loop sums them in, so each element takes the same terms in the same order) regroup each operand as a stack of matrices, A as batch × free_a × summed and B as batch × summed × free_b, multiplied batch by batch through the packed kernel, batches in parallel; the result is laid out in the output's order; a dot product per batch by dot_product. Nothing for a label repeated within an operand or summed by one operand alone (the general loop takes those). | |
| std::string | pair_keeps (const std::vector< std::string > &labels, std::size_t i, std::size_t j, const std::string &output) |
The labels of labels[i] and labels[j] that the output or another operand carries, in their order of appearance, once each: what their contraction keeps. | |
| double | extent_product (const std::string &labels, const std::map< char, std::size_t > &sizes) |
The product of the extents of labels' distinct labels, as a double (a cost: it may pass size_t). | |
| std::vector< std::string > | after_step (const std::vector< std::string > &labels, contraction_step step, const std::string &output) |
| The list after a contraction step: the two operands' labels removed, their result's appended. | |
| void | cheapest_path (const std::vector< std::string > &labels, const std::string &output, const std::map< char, std::size_t > &sizes, std::vector< contraction_step > &path, double &cost, double &largest) |
The cheapest order of contraction from labels, and its cost and largest intermediate: every order tried, the first found kept among equals. | |
| std::vector< contraction_step > | contraction_path (std::vector< std::string > labels, const std::string &output, const std::map< char, std::size_t > &sizes) |
The order in which to contract operands labelled labels into output, as opt_einsum chooses it: the cheapest in multiply-adds (the product of the extents a pair's labels span), the smaller largest intermediate among equals, over every order up to exhaustive_path_operands operands; past it, the cheapest pair at each step among those sharing a label (an outer product, cheap in itself, would leave its labels to every later step), the smaller result, then the first pair, among equals. It depends on the shapes alone, so a result is the same on any number of threads. | |
| std::size_t | ellipsis_in (const std::string &subscripts) |
| Evaluates a parsed einsum over its operands. | |
| template<typename T > | |
| tensor< T > | einsum_impl (const std::string &spec_in, const std::vector< const tensor< T > * > &operands) |
| template<typename T > | |
| tensor< T > | einsum_ellipsis (const std::string &spec, const std::vector< const tensor< T > * > &operands) |
einsum_impl of a spec with ellipses, as NumPy reads them: each ... stands for an operand's axes its letters leave, aligned from the right across the operands and broadcast (an extent of 1 repeated), named by letters the spec does not use. An operand's axis of extent 1 under the ellipsis where another's is longer is dropped with its letter, so that operand repeats along it. The output: its ... the broadcast axes (an output without one, where an operand has such axes, is refused, as NumPy refuses it); with no ->, the broadcast axes then the letters that appear once, in order. | |
| std::string | shape_text (const std::vector< std::size_t > &shape) |
A shape as NumPy prints it, (2, 3), for a message. | |
| std::vector< std::size_t > | broadcast_strides (const std::vector< std::size_t > &shape, std::size_t rank) |
The row-major strides of shape, with a stride of 0 on each axis of extent 1 (it broadcasts), aligned on the last axis of a result of rank rank (missing leading axes stride 0 too). | |
| axis_walk | walk_along (const std::vector< std::size_t > &shape, std::size_t axis) |
axis_walk of shape along axis; outer and inner 0 where the result is empty. | |
| template<std::floating_point R> | |
| R | numpy_muladd (R x, R y, R z, bool fused) |
x · y + z: one fused multiply-add where fused, else the product rounded first (statements apart, so no compiler contracts them). | |
| template<split_complex C> | |
| C | numpy_multiply (const C &a, const C &b, bool fused) |
a · b of complex numbers as NumPy's array loop computes it: ai·bi and ai·br rounded, then ar·br − ai·bi and ar·bi + ai·br, each by a fused multiply-add where fused (NumPy's vectors on arm64 and on x86-64 with FMA), else each product rounded first. No recovery of an infinity from a NaN, as std::complex's (C Annex G) makes: 1j · (inf + inf·j) is NaN, (1e308 + 1e308j) · (1e308 − 1e308j) inf + inf·j when fused. | |
| template<split_complex C> | |
| C | numpy_divide (const C &a, const C &b, bool fused) |
a / b of complex numbers as NumPy's array loop computes it: Smith's division by the larger part of b (an exact zero divided part by part, to an infinity or a NaN), each x + y·z of it one fused multiply-add where fused (NumPy's build contracts them on arm64, not on x86-64). | |
| bool | numpy_fuses_products () noexcept |
| Whether NumPy's complex product fuses its multiply-adds on this processor: on arm64, and on x86-64 where its AVX2 loop runs (FMA with it); elsewhere its loop rounds each product. | |
| bool | numpy_fuses_scalar_products () noexcept |
| Whether NumPy's complex product fuses its multiply-adds in its scalar loop (where numpy_loop_is_scalar holds), as its PyPI wheels are built: that loop sits beside the vector one, and clang (macOS) and MSVC (Windows) contract it where the vector loop has FMA; GCC (Linux, ISO C) does not. | |
| bool | numpy_fuses_quotients () noexcept |
| Whether NumPy's complex quotient fuses its multiply-adds, as its PyPI wheels are built: on macOS arm64, whose clang contracts them (FMA is arm64's baseline); not on Linux, whose GCC does not, nor on x86-64, whose baseline has no FMA. | |
| template<typename T > | |
| bool | numpy_loop_is_scalar (const tensor< T > &a, const tensor< T > &b) noexcept |
Whether NumPy multiplies a by b in its scalar loop rather than its vector one: when the result is one element and an operand of rank above 0 is broadcast to a higher rank, NumPy's iterator hands its loop an output step of 0, which the vector loop refuses. | |
| template<std::floating_point T> | |
| T | divided_by_count (T value, std::size_t n) |
A mean: value divided by the count n, in the wider of T and double and back to T, the count exact (to 2^53), as NumPy divides by its integer count; a float's count rounded to a float past 2^24 (2^24 + 1 ones averaged to 1, NumPy's 0.99999994). | |
| template<std::floating_point T, typename Order > | |
| T | settled_zero (T value, const T *first, std::size_t count, std::size_t stride, const Order &larger) |
The extreme value of count elements stride apart from first, its zero's sign settled as IEEE 754-2019's maximum and minimum settle it, and NumPy's on arm64: no zero is larger than the other, so the first met stayed; max gives +0 where one is among them, min −0, whatever their order. A NaN, or any other value, as it is. | |
| std::vector< std::size_t > | without_axis (std::vector< std::size_t > shape, std::size_t axis) |
The shape with axis axis removed. | |
| template<typename T , typename Reduce > | |
| tensor< T > | reduce_along (const tensor< T > &t, std::size_t axis, const Reduce &reduce) |
Reduces each line of t along axis by reduce, which takes the line's first element, its length and its stride. | |
| template<typename T > | |
| T | square (const T &x) |
x · x. | |
| template<typename U > | |
| std::complex< U > | square (const std::complex< U > &x) |
| NumPy's complex square, (re·re − im·im) + (re·im + im·re)i, without the infinity handling of a general product (C Annex G), which costs a call per element. | |
| template<typename T > | |
| bool | is_nan (const T &) |
| A NaN? Never for a scalar without one. | |
| bool | is_nan (double x) |
| A NaN? | |
| template<typename T , typename Order > | |
| T | contiguous_extreme (const T *first, std::size_t count, const Order &larger) |
The largest (larger is std::greater) or smallest of count contiguous elements, or a NaN if one is among them, as NumPy's: the comparisons run without a branch over native vectors (or eight scalar lanes), each lane its own running best, and the NaN is tracked apart. | |
| template<typename T , typename Order > | |
| T | extreme (const T *first, std::size_t count, std::size_t stride, const Order &larger) |
The largest (larger is std::greater) or smallest element of a line. | |
| template<typename T , typename Order > | |
| tensor< T > | extremes_along (const tensor< T > &t, std::size_t axis, const Order &larger) |
The largest (larger is std::greater) or smallest elements along axis. Along an inner axis each line is contiguous (extreme); along an outer one the slabs are compared element by element, so the loop runs over contiguous memory rather than down a column. | |
| std::size_t | checked_mul (std::size_t a, std::size_t b) |
a * b, throwing numerical::dimension_error if the product overflows std::size_t. Guards dimension/shape products (element counts, allocation sizes) against silent wraparound on absurd inputs. The a != 0 guard also keeps the max / a division well-defined. | |
| template<typename It > | |
| std::size_t | element_count (It first, It last) |
The elements a shape's extents [first, last) hold: 0 where one of them is 0, whatever the others, else their product (checked_mul). The extents a tensor keeps, past one it holds 0 of, may multiply past std::size_t: a result of them is empty, not of a wrapped size. | |
Variables | |
| constexpr std::size_t | 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 | 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 | lu_small_block {32} |
| The columns of one panel of the LU of an order up to lu_small_order. | |
| constexpr std::size_t | 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 | 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 | 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 | 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 | solve_leaf {16} |
| The columns of one leaf of the recursive LU panel. | |
| constexpr std::size_t | 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 | 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 | svd_band {32} |
| The rows of an SVD's input one thread scales and copies at a time. | |
| constexpr std::size_t | zeroed_rows {32} |
| The rows of a Cholesky factor whose upper part one thread zeroes at a time. | |
| constexpr std::size_t | cholesky_block {64} |
| The columns of one panel of the blocked Cholesky factorisation (as lu_block). | |
| constexpr std::size_t | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | plain_squares_high {0x1p900} |
| The largest plain sum of squares plain_or_at_scale keeps: no term of it overflowed. | |
| constexpr std::size_t | finite_chunk {std::size_t {1} << 15U} |
| The elements all_finite_shared gives a thread at a time. | |
| constexpr std::size_t | qr_block {32} |
| The reflectors of one panel of the blocked QR. | |
| constexpr std::size_t | 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 | 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 | 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 | tridiagonal_block {32} |
| The reflections of one panel of the blocked tridiagonal reduction. | |
| constexpr std::size_t | 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 | 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 | divide_leaf {32} |
| The order up to which tridiagonal_divide_and_conquer diagonalises by implicit QL. | |
| constexpr std::size_t | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | 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 | bidiagonal_leaf {32} |
| The order up to which bidiagonal_divide decomposes by implicit-shift QR. | |
| constexpr double | 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. | |
| constexpr double | schur_ulp {std::numeric_limits<double>::epsilon()} |
| The relative spacing LAPACK's ‘dlamch('P’)` names: the roundoff a QR sweep judges a subdiagonal by. | |
| constexpr std::size_t | schur_steps_per_order {30} |
The iterations one eigenvalue (or pair) may take per unit of the active order, as dlahqr's 30. | |
| constexpr double | balance_limit {0x1p256} |
| The largest power of two a balancing scales a row or a column by, either way: the matrix enters at a largest magnitude near 1, so no entry overflows, and the scaling stops where LAPACK's would be one factor among dozens. | |
| constexpr double | back_substitution_limit {0x1p400} |
| The modulus an eigenvector's components are kept under during the back-substitution: the sums stay finite for any order a matrix can have, and a component that would grow past it scales the vector down. | |
| constexpr std::size_t | hessenberg_block {32} |
| The columns of a panel of hessenberg_blocked. | |
| constexpr std::size_t | hessenberg_crossover {128} |
| The active order from which hessenberg reduces by panels. | |
| constexpr lapack_index | multishift_crossover {75} |
LAPACK's nmin: the active order below which dlaqr0's iterations skip a sweep after enough deflation. | |
| constexpr lapack_index | window_growth {5} |
The iterations without deflation after which the window grows (kexnw). | |
| constexpr lapack_index | exceptional_after {6} |
The iterations without deflation after which the shifts are exceptional (kexsh). | |
| constexpr lapack_index | nibble {14} |
The share of the window, in percent, whose deflation makes the sweep unnecessary (nibble). | |
| constexpr bool | neon_compiled {false} |
| Whether the NEON path's tile is compiled in: on every arm64 target GCC or Clang builds. | |
| constexpr std::size_t | gemm_mr {4} |
| The rows of C one register tile covers. | |
| template<typename T > | |
| constexpr std::size_t | gemm_nr {sizeof(T) >= 16 ? std::size_t {6} : std::size_t {96} / sizeof(T)} |
The columns of C one register tile covers: 96 bytes of elements (12 double, 24 float), which with gemm_mr rows fills the vector registers of both arm64 and x86-64 with AVX. | |
| constexpr std::size_t | gemm_kc {256} |
The depth of one panel: a B sliver (gemm_kc × gemm_nr) stays in the L1 cache. | |
| constexpr std::size_t | gemm_mc {64} |
The rows of one task: an A block (gemm_mc × gemm_kc) stays in the L2 cache. | |
| constexpr std::size_t | gemm_slivers_per_task {4} |
| The B slivers of one task. | |
| constexpr std::size_t | tasks_per_thread {4} |
| The tasks a product's compute pass hands each thread at least, while its rows allow: a slower core then holds back no one. | |
| constexpr std::size_t | gemm_kernel_speedup {16} |
| The multiply-adds the register tile does in the time of one in work_per_thread's dependent chain (about 20 G/s on one arm64 or x86-64 AVX2 core): a task weighs its work by it, so a product of 0.1 ms does not start eight threads. | |
| template<typename T > | |
| constexpr std::size_t | gemm_nc {gemm_nr<T> * 128} |
| The columns of B packed at once, a whole number of slivers. | |
| constexpr std::size_t | gemm_packed_threshold {std::size_t {1} << 13U} |
| Below this many multiply-adds the plain loop (gemm_rows, by blocks of rows held in registers) runs instead of the packed kernel, both summing in the same order. Against the packed kernel on 400 random shapes from 512 multiply-adds, the plain loop took 0.25 to 0.41 times as long below 2^13 on arm64 (M1 Pro, geometric means by octave, the slowest shape 0.91); above, it lost up to 2.5 times on shapes the AMX suits. On x86-64 the plain loop has the baseline's vectors and the packed one the AVX2 or AVX-512 units: there it took 0.43 to 0.55 times as long below 2^11, 0.86 above with shapes up to 2 times slower (six AVX2 cores, GCC 13). | |
| constexpr std::size_t | gemm_min_depth {8} |
| Below this depth the plain loop runs too: each panel reads and writes its tiles of C once, which a shallow product does not amortise (1000×1000×1 ran six times slower packed, ×4 still lost 10 %). | |
| constexpr std::size_t | narrow_least_depth {2} |
| The least depth at which rows of C wider than one block (narrow_columns) run by blocks in gemm_rows; one block or narrower takes any depth past 1. Against the plain loop of unknown width, on 500 random shapes below the packed kernel, the blocks took 0.11 to 0.63 times as long from a depth of 2 on arm64 (M1 Pro, geometric means by width); on x86-64, where the loop has the baseline's vectors, 0.35 to 0.90 from a depth of 8 but 1.08 to 1.31 below it past one block (six AVX2 cores, GCC 13). | |
| constexpr std::size_t | narrow_columns {8} |
| The columns of C one block of gemm_rows keeps in registers over the whole depth. | |
| constexpr std::size_t | narrow_group {4} |
| The rows of C one block of gemm_rows takes together: each element a chain of multiply-adds that waits on its last, four rows give the units four times the chains (M1 Pro, order 8 to 24: 1.6 to 1.9 times as fast as one row; two rows on six x86-64 AVX2 cores, the fastest there, were within 7 %). | |
| constexpr std::size_t | packed_line {64} |
| The cache line a packed operand starts. | |
| template<typename Kernel > | |
| constexpr bool | prefetches_next_row_of_c |
Whether Kernel has the C tile of the next row of tiles prefetched while a tile computes (Kernel::prefetch_rows): the tiles within a row walk one page's lines, which the processor's own prefetch follows, while the next row lands on other pages. | |
| constexpr std::size_t | split_columns {512} |
The columns of B and C, and the depth, of one real product of gemm_split: its real operands stay a few megabytes whatever the shapes. | |
| constexpr std::size_t | split_depth {256} |
The depth of one real product of gemm_split (split_columns). | |
| constexpr std::size_t | split_kept {std::size_t {1} << 18} |
The most elements an operand of gemm_split keeps from one product to the next on a thread: 2 MiB of double, every operand of a product of order 512; a larger one is let go after its product. | |
| constexpr std::size_t | split_parts_rows {64} |
The most rows of A for which gemm_split runs four real products of a block's parts on the AMX path (split_by_parts). M1 Pro, all threads, A of 16 to 64 rows by B of 256 to 1024: complex float 1.1 to 2.4 times faster, complex double 0.97 to 1.6; 96 to 128 rows by 512, complex double 1.2 to 1.4 times slower, where one wide product suited the AMX. | |
| constexpr std::size_t | plain_loop_least_width {64} |
| The fewest columns of a product plain_loop_pays sends to the plain loop on the x86-64 tiles: there an order-24 SVD, whose two products of 24 × 24 × 24 it took, ran 9 % slower (GCC 13, the same bits), though such a product alone ran faster in the plain loop. On the NEON tile the small square products gain most (GCC 14, one core: a product of order 24 to 64 −32 to −50 %, a QR −5 to −21 %). | |
| constexpr std::size_t | dot_chunk {std::size_t {1} << 14U} |
The elements of one partial sum of chunked_dot: fixed, so the sum does not depend on the threads. | |
| template<typename T > | |
| constexpr std::size_t | multiply_add_weight {1} |
The real multiply-adds one multiply-add of T weighs: 1 for a real number. | |
| constexpr std::size_t | row_group {64} |
| The rows a thread takes of a matrix times a vector (matrix_times_vector). | |
| constexpr std::size_t | no_core {std::numeric_limits<std::size_t>::max()} |
| No core: a worker free to run anywhere, or a caller on a processor the pool does not know. | |
| constexpr std::chrono::microseconds | spin_time {500} |
| How long a thread polls before it sleeps: kernels run in quick succession (a matrix product packs, then computes, panel after panel), and a thread woken from sleep starts tens of microseconds late, its wake-up paid by the caller of the next run, about 10 µs a sleeping worker on six x86-64 cores. A factorisation's steps last about 100 µs at order 256, so at 100 µs the workers one step left out were asleep for the next (an LU of order 256: six x86-64 cores 560 → 500 µs at 500; eight arm64 cores 358 → 324). | |
| constexpr std::size_t | work_per_thread {std::size_t {1} << 16} |
| The work one thread must have to pay for starting it: about 60 µs of multiply-adds, six times the cost of a start (measured on arm64: 16 indices of 4096 ran 80 µs sequentially, 92 µs on eight threads; of 16384, 324 µs and 128 µs). | |
| constexpr std::size_t | elementwise_chunk {work_per_thread / 4} |
| The elements one task of an element-wise operation takes, a quarter of a thread's work: a tensor of a million shares out (M1 Pro, 2000 × 2000 doubles added: 3.1 → 1.4 ms), one of a few hundred thousand stays within its time (a sixteenth: 300 × 300 plus a row 13 % slower, each task setting up its odometer). | |
| constexpr std::size_t | vector_speedup {4} |
| The multiply-adds a vectorised inner loop does in the time of one in work_per_thread's dependent chain: a product weighs its work by it (at order 64, a matrix product on eight threads took 95 µs, on one 35 µs). | |
| constexpr std::size_t | zero_chunk {std::size_t {1} << 12} |
The elements one index of parallel_zero writes: eight 4 KiB pages of double. | |
| constexpr std::size_t | zero_work {1} |
The work of zeroing one element, in multiply-adds: a second thread from 1 MiB of double. Below, the allocator mostly hands back memory already touched, and zeroing 512 KiB took 4.5 µs on one arm64 thread, 9 µs split across eight. | |
| constexpr std::size_t | transpose_tile {32} |
The side of the tiles a transposing copy goes by (an einsum permutation of two axes, the eigenvectors' columns): a tile read and a tile written fit the L1 cache, where element by element each write, a row apart, missed it (order 1024 on arm64: 4.7 → 1.5 ms; the eigenvectors of order 1024, a tenth of an eigh on six x86-64 cores). | |
| constexpr std::size_t | symmetric_row_group {4} |
| The rows of a lower triangle symmetric_rows takes at once. | |
| constexpr std::size_t | product_rows {8} |
The rows of y subtract_rows takes at once. | |
| constexpr std::size_t | pairwise_columns {512} |
| The columns pairwise_sum_rows sums at once: its partial sums of eight rows, 32 KiB of doubles, stay in the L1 cache. | |
| constexpr std::size_t | vector_chunk {128} |
| The columns a thread takes of a vector times a matrix: each of its rows added to them in order. | |
| constexpr std::size_t | transposed_columns {8} |
| The columns of B one task of vector_times_transposed sums together: their sums' chains interleaved, each waiting on its own last addition. | |
| constexpr std::size_t | exhaustive_path_operands {4} |
| The operands at most whose contraction order contraction_path searches exhaustively; more take the cheapest pair at each step. | |
| using numerical::detail::chunk_runner = typedef void (*)(const void* body, std::size_t begin, std::size_t end) |
A chunk of a parallel_for, type-erased: runs body(i) for i in [begin, end).
Definition at line 234 of file parallel.hpp.
| using numerical::detail::lapack_index = typedef std::ptrdiff_t |
The type std::abs gives for a T: the real type of a complex one, T itself otherwise.
Definition at line 268 of file decomposition.hpp.
| using numerical::detail::rounding_is = typedef std::integral_constant<term_rounding, Mode> |
The rounding mode Mode as a type, which picks rounded_term's overload.
A task the caller of a run does before it takes chunks (parallel_for_beside), type-erased.
Definition at line 237 of file parallel.hpp.
| using numerical::detail::symmetric_row_kernel = typedef 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.
Definition at line 2729 of file decomposition.hpp.
| using numerical::detail::symmetric_rows_kernel = typedef 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.
Definition at line 2732 of file decomposition.hpp.
| using numerical::detail::value_range_t = typedef 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).
Definition at line 295 of file decomposition.hpp.
|
strong |
The tile a product's left operand is packed for (numerical::detail::pack_left): plain where it is not packed.
| Enumerator | |
|---|---|
| plain | gemm_rows, or a scalar whose products are not shared: A is read in place. |
| vector | The portable path's vector tile. |
| neon | The NEON tile. |
| amx | The AMX tile. |
| avx512 | The AVX-512 tile. |
| avx2 | The AVX2 tile. |
How numerical::detail::gemm_rows rounds each term sum + a·b, as the tile of the path in use does: as the compiler rounds it (the vector tile's way), one fused multiply-add (arm64's NEON and AMX tiles), or one fused multiply-add of AVX2's unit (the x86-64 tiles), in code built for any x86-64 processor.
| Enumerator | |
|---|---|
| compiler | As the compiler rounds |
| fused |
|
| x86_fused |
|
| 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.
Definition at line 2692 of file decomposition.hpp.
| 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).
Definition at line 2707 of file decomposition.hpp.
| void numerical::detail::add_one_product | ( | C & | y, |
| const C & | f, | ||
| const C & | x | ||
| ) |
add_one_product of complex numbers: complex_scaled_add, as add_product computes them.
y += f · x of one element, as add_product computes each of its elements.
| void numerical::detail::add_product | ( | C * | y, |
| const C & | f, | ||
| const C * | x, | ||
| std::size_t | count | ||
| ) |
add_product of complex numbers: vectors of their parts, then complex_scaled_add, every element the same arithmetic.
|
inline |
y[i] += f · x[i] for double, as subtract_product.
|
inline |
add_product for complex double (avx2_scaled_add_c64 on the x86-64 paths).
|
inline |
y[i] += f · x[i] for complex float: on the x86-64 paths by AVX2 vectors of its arithmetic (avx2_scaled_add_c32), the same bits as add_product's template, which runs elsewhere.
y[i] += f · x[i], as subtract_product.
The conjugate transpose of a.
Definition at line 5320 of file decomposition.hpp.
|
inlineprivate |
The list after a contraction step: the two operands' labels removed, their result's appended.
Definition at line 734 of file tensor.hpp.
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.
Definition at line 505 of file decomposition.hpp.
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).
Definition at line 560 of file decomposition.hpp.
|
inlinenoexcept |
|
private |
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.
Definition at line 2027 of file decomposition.hpp.
| 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.
A run of chained rotations goes as LAPACK's dlasr: the row each rotation hands the next stays in a buffer the size of the chunk, so every row of the run is read and written once rather than twice per rotation.
Definition at line 3408 of file decomposition.hpp.
|
private |
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).
Definition at line 2070 of file decomposition.hpp.
|
inline |
The indices of diagonal by ascending value.
Definition at line 4344 of file decomposition.hpp.
Balances the row-major n × n a in place, as LAPACK's ‘dgebal('B’)`: rows with no off-diagonal entry in the active block go to its bottom and columns with none to its top, each isolating an eigenvalue on the diagonal; then each row and column of the remaining block is scaled by a power of two until their 2-norms are within a factor of 2, which leaves every entry exact.
|
inlineprivate |
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.
The rows split at a middle row k: above it an k × (k + 1) bidiagonal, below it an (n − k − 1) × (…) one, each decomposed the same way (by implicit-shift QR up to bidiagonal_leaf); row k's two entries, in their right bases, with the two null vectors rotated into one, make an arrow, which merge_arrow decomposes.
The tree is laid out first (plan_bidiagonal). The blocks at the depth that gives each thread one, and the leaves above it, are independent of one another and go to the threads together, each decomposed on one thread, a block right after its halves while they are still in its caches; the joins above them then come a level at a time, those of one level together, a level of one join sharing it among the threads. Every block takes the same inputs either way, so the result is the same to the bit; on one thread the order is the recursion's.
Definition at line 5918 of file decomposition.hpp.
|
private |
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.
| convergence_error | If an iteration did not converge within its bound. |
Definition at line 5966 of file decomposition.hpp.
|
private |
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.
| [in,out] | d | The diagonal; on return, the singular values, descending and non-negative. |
| [in,out] | e | The superdiagonal in e[0…n-2], e[n-1] = 0; destroyed. |
| [in,out] | ut | Row-major n × n; on return, row j is the left vector of d[j]. Empty, with vt: the singular values alone. |
| [in,out] | vt | Row-major n × n; on return, row j is the right vector of d[j]. |
| [in] | n | The order. |
| [in] | max_steps | The bound of QR steps in all. |
| convergence_error | If the iteration has not converged after max_steps steps. |
Definition at line 5174 of file decomposition.hpp.
|
private |
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.
| [in,out] | b | The matrix; destroyed. |
| [in] | n | Its order. |
| [out] | d | The diagonal (size n). |
| [out] | e | The superdiagonal in e[0…n-2], and e[n-1] = 0 (size n). |
| [out] | ut | U_bᵀ, row-major: its rows are the columns of U_b (left as it is without form). |
| [out] | vt | V_bᵀ, row-major: its rows are the columns of V_b (likewise). |
| [in] | form | Whether to form ut and vt; the bidiagonal alone gives the singular values. |
Definition at line 4871 of file decomposition.hpp.
|
private |
bidiagonalize_unblocked, for a complex or any other T the panels do not take.
Definition at line 5149 of file decomposition.hpp.
|
private |
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).
A complex entry left on the bidiagonal is turned real by a unit phase on its row or column, as LAPACK does, so the iteration that follows needs only real rotations.
| [in,out] | b | The matrix; destroyed. |
| [in] | n | Its order. |
| [out] | d | The diagonal (size n). |
| [out] | e | The superdiagonal in e[0…n-2], and e[n-1] = 0 (size n). |
| [out] | ut | U_bᵀ, row-major: its rows are the columns of U_b (left as it is without form). |
| [out] | vt | V_bᵀ, row-major: its rows are the columns of V_b (likewise). |
| [in] | form | Whether to form ut and vt. |
Definition at line 5078 of file decomposition.hpp.
|
inline |
|
inline |
The row-major strides of shape, with a stride of 0 on each axis of extent 1 (it broadcasts), aligned on the last axis of a result of rank rank (missing leading axes stride 0 too).
Definition at line 48 of file tensor_ops.hpp.
|
private |
The product of the row-major rows × inner a and inner × cols b, shaped out_shape: the canonical two-operand contraction of einsum_impl (a matrix times a vector, a vector times a matrix, a matrix product, a dot product), on any storage that holds the operands contiguously.
Definition at line 399 of file tensor.hpp.
| void numerical::detail::chase_bulges | ( | const schur_target< T > & | target, |
| lapack_index | ktop, | ||
| lapack_index | kbot, | ||
| std::complex< double > * | s, | ||
| lapack_index | count, | ||
| bool | accumulate | ||
| ) |
One small-bulge multishift QR sweep over the active block ktop … kbot of target's Hessenberg matrix with the count shifts at s, as LAPACK's dlaqr5 and zlaqr5: chains of 3 × 3 bulges, two shifts each, chased down the diagonal by slabs of 4 · count / 2 columns; with accumulate, each slab's reflections are accumulated in U and applied to the rest of the matrix and to Z through the matrix product kernel, else applied to them directly (LAPACK's kacc22, which accumulates from 14 shifts). A subdiagonal entry that becomes negligible on the way is set to zero (vigilant deflation).
|
inlineprivate |
The cheapest order of contraction from labels, and its cost and largest intermediate: every order tried, the first found kept among equals.
Definition at line 750 of file tensor.hpp.
|
inline |
a * b, throwing numerical::dimension_error if the product overflows std::size_t. Guards dimension/shape products (element counts, allocation sizes) against silent wraparound on absurd inputs. The a != 0 guard also keeps the max / a division well-defined.
Definition at line 50 of file vector.hpp.
|
constexprnoexcept |
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).
Definition at line 223 of file decomposition.hpp.
Σ a[i] · b[i] for a plain dot product: chunks of dot_chunk elements in parallel, each by dot_product (several vector accumulators, on the widest x86-64 units for double), their sums added in the chunks' order. The same on any number of threads; in the last bits, not the ascending sum, which waits on each addition before the next (ten times slower at 100 000 elements).
| 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.
Definition at line 4360 of file decomposition.hpp.
| void numerical::detail::complex_row_dots | ( | const C * | m, |
| std::size_t | ld, | ||
| std::size_t | rows, | ||
| const C * | x, | ||
| std::size_t | count, | ||
| C * | out | ||
| ) |
row_dots of complex numbers: each row by complex_dot_product, as chunked_dot sums a row.
| void numerical::detail::complex_scaled_add | ( | C & | y, |
| const C & | f, | ||
| const C & | x | ||
| ) |
y += f · x of one complex number, its part products rounded apart as libc++'s operator computes them (add_scaled_vectors's arithmetic, for the elements past its vectors).
|
inline |
The complex Schur form of the upper Hessenberg t's active block, as LAPACK's zlahqr: the shifted QR algorithm with Wilkinson's shift, the deflation criterion of Ahues and Tisseur, and exceptional shifts after 10 and 20 iterations without deflation. The eigenvalues are left in values.
| [in] | t | The Hessenberg matrix, its active block, and the Schur vectors to update. |
| [in,out] | values | The eigenvalues (the isolated ones already in place). |
| [in] | steps | The iterations allowed between two deflations. |
| convergence_error | If an eigenvalue has not converged within steps. |
The conjugate of a complex scalar.
Definition at line 261 of file decomposition.hpp.
The conjugate of a real (or general) scalar — the value itself.
Definition at line 254 of file decomposition.hpp.
Σ conj(x[i]) · y[i] for a real T, where the conjugate is the value: dot_product.
Definition at line 2105 of file decomposition.hpp.
Σ conj(x[i]) · y[i] of complex numbers (conjugate_dot_product).
Definition at line 2114 of file decomposition.hpp.
Σ conj(x[i]) · y[i], in order, where no native vector holds the parts.
| T numerical::detail::contiguous_extreme | ( | const T * | first, |
| std::size_t | count, | ||
| const Order & | larger | ||
| ) |
The largest (larger is std::greater) or smallest of count contiguous elements, or a NaN if one is among them, as NumPy's: the comparisons run without a branch over native vectors (or eight scalar lanes), each lane its own running best, and the NaN is tracked apart.
Definition at line 380 of file tensor_ops.hpp.
|
private |
A two-operand contraction as matrix products: the labels both operands and the output carry (batch), those one operand and the output carry (free), and those both operands carry alone (summed, in ascending order, the order the general loop sums them in, so each element takes the same terms in the same order) regroup each operand as a stack of matrices, A as batch × free_a × summed and B as batch × summed × free_b, multiplied batch by batch through the packed kernel, batches in parallel; the result is laid out in the output's order; a dot product per batch by dot_product. Nothing for a label repeated within an operand or summed by one operand alone (the general loop takes those).
Definition at line 610 of file tensor.hpp.
|
inlineprivate |
The order in which to contract operands labelled labels into output, as opt_einsum chooses it: the cheapest in multiply-adds (the product of the extents a pair's labels span), the smaller largest intermediate among equals, over every order up to exhaustive_path_operands operands; past it, the cheapest pair at each step among those sharing a label (an outer product, cheap in itself, would leave its labels to every later step), the smaller result, then the first pair, among equals. It depends on the shapes alone, so a result is the same on any number of threads.
Definition at line 794 of file tensor.hpp.
| 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.
| [in] | block | The block's first element; its rows n apart. |
| [in] | n | The row stride. |
| [in] | m | The block's order. |
| [in] | v | The vector (m elements). |
| [out] | w | B v corrected (m elements). |
| [in,out] | scratch | hermitian_product's partials, kept by the caller from one step to the next. |
| [in] | corrections | The rows c₀ … c_{terms−1} (ldc apart), over the block's rows. |
| [in] | ldc | The distance between those rows. |
| [in] | terms | Their count, even (two a step: its v and its w). |
| [out] | factors | The coefficients (terms). |
| [in,out] | next | The panel's next column over the block's rows, which takes the terms of these corrections with next_factors in the same pass; none past a panel's last step. |
| [in] | next_factors | Its coefficients of the rows (terms). |
Definition at line 3018 of file decomposition.hpp.
|
inline |
|
inlinenoexcept |
This process's identifier: a pool started before a fork belongs to the parent, whose threads the child does not have.
Definition at line 49 of file parallel.hpp.
|
inline |
Cyclic Jacobi on the full symmetric a, the rotations accumulated into v.
| [in,out] | a | The matrix; on return, diagonal to rounding. |
| [in,out] | v | The rotations (start from the identity). |
| [in] | n | The order. |
| [in] | frobenius | The Frobenius norm of a, the scale of the stopping threshold. |
| [in] | max_sweeps | The bound of sweeps. |
| convergence_error | If the off-diagonal is not negligible after max_sweeps sweeps. |
Definition at line 4322 of file decomposition.hpp.
|
inline |
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.
The tree is laid out first (plan_divide). The blocks at the depth that gives each thread one, and the leaves above it, are independent of one another and go to the threads together, each solved on one thread with its own buffers, a block right after its halves while they are still in its caches; the joins above them then come a level at a time, those of one level together, a level of one join sharing it among the threads. One after another each block had run alone (six AVX2 cores, an eigh of order 256: the leaves and the joins below the last, 1.5 ms of the 4 ms the eigenvectors took, no faster than on one). Every block takes the same inputs either way, so the result is the same to the bit; on one thread the order is the recursion's.
Definition at line 4144 of file decomposition.hpp.
| scaled_quotient< S > numerical::detail::divided | ( | S | r, |
| P | pivot, | ||
| double | floor | ||
| ) |
r / pivot, a pivot below floor taken as floor, r scaled first where the quotient would pass back_substitution_limit.
A mean: value divided by the count n, in the wider of T and double and back to T, the count exact (to 2^53), as NumPy divides by its integer count; a float's count rounded to a float past 2^24 (2^24 + 1 ones averaged to 1, NumPy's 0.99999994).
divided_by_count of any other scalar: divided by the count as a T.
divided_by_count of a complex number of float or double, as NumPy's mean divides it: by the complex count n + 0j, in complex double, by its division (numpy_divide), so a NaN part makes both NaN and each part is multiplied by the reciprocal of the count.
| dimension_error | If the count is 0: an integer has no NaN for the mean of nothing. |
Definition at line 213 of file tensor_ops.hpp.
|
inline |
Σ x[i] · y[i] for double, as subtract_product chooses its units (in another order of summation there, as fixed).
|
inline |
dot_product for std::complex<double>, through complex_dot_product.
|
inline |
dot_product for std::complex<float>, through complex_dot_product.
| deflation_count numerical::detail::early_deflation | ( | const schur_target< T > & | target, |
| lapack_index | ktop, | ||
| lapack_index | kbot, | ||
| lapack_index | nw, | ||
| std::complex< double > * | w, | ||
| bool | recursive | ||
| ) |
Aggressive early deflation, as LAPACK's dlaqr3 and zlaqr3: the trailing nw × nw window of the active block ktop … kbot is brought to Schur form; its eigenvalues whose coupling to the rest (the spike, s times the first row of the window's Schur vectors) is negligible deflate; the window is reduced back to Hessenberg form above them, and its transformation carried to the rest of H and to Z through the matrix product kernel. Every eigenvalue of the window is written to w: the deflated ones final, the others the next shifts. recursive: the window itself by multishift_schur past its crossover.
The magnitudes of the eigenvalues of the real symmetric or complex Hermitian a, its lower triangle read: its singular values, unordered.
Definition at line 6724 of file decomposition.hpp.
|
private |
einsum_impl of a spec with ellipses, as NumPy reads them: each ... stands for an operand's axes its letters leave, aligned from the right across the operands and broadcast (an extent of 1 repeated), named by letters the spec does not use. An operand's axis of extent 1 under the ellipsis where another's is longer is dropped with its letter, so that operand repeats along it. The output: its ... the broadcast axes (an output without one, where an operand has such axes, is refused, as NumPy refuses it); with no ->, the broadcast axes then the letters that appear once, in order.
Definition at line 874 of file tensor.hpp.
|
private |
Definition at line 986 of file tensor.hpp.
The elements a shape's extents [first, last) hold: 0 where one of them is 0, whatever the others, else their product (checked_mul). The extents a tensor keeps, past one it holds 0 of, may multiply past std::size_t: a result of them is empty, not of a wrapped size.
| dimension_error | If none is 0 and their product overflows. |
Definition at line 66 of file vector.hpp.
|
inlineprivate |
Evaluates a parsed einsum over its operands.
After parsing and the consistency checks, the contraction runs over a single integer index space — the output (free) labels first, then the summed labels. Each operand's flat offset is carried incrementally by a mixed-radix odometer using a per-label stride coefficient computed once, so the hot loop allocates nothing and performs no map lookups; its innermost reduction is a strided accumulation the compiler can auto-vectorise. (Pairwise contraction-path optimisation, as in opt_einsum, and a BLAS reduction are deliberately not included — they would trade the library's generic simplicity and zero dependencies for raw throughput.)
The position of the ellipsis in subscripts, or npos where there is none: ‘’.'only as one.... \throws einsum_error If a'.'` stands otherwise.
Definition at line 851 of file tensor.hpp.
|
inline |
|
inline |
|
inline |
Exchanges the adjacent diagonal blocks of the real quasi-triangular t at j1, of orders n1 then n2, by an orthogonal similarity carried to the rows of qt, as LAPACK's dlaexc.
10 ε ‖D‖ is refused and nothing changes. | void numerical::detail::exchange_diagonal | ( | one_based< T > | t, |
| one_based< T > | qt, | ||
| lapack_index | n, | ||
| lapack_index | k | ||
| ) |
|
inlineprivate |
The product of the extents of labels' distinct labels, as a double (a cost: it may pass size_t).
Definition at line 719 of file tensor.hpp.
| T numerical::detail::extreme | ( | const T * | first, |
| std::size_t | count, | ||
| std::size_t | stride, | ||
| const Order & | larger | ||
| ) |
The largest (larger is std::greater) or smallest element of a line.
Definition at line 437 of file tensor_ops.hpp.
| tensor< T > numerical::detail::extremes_along | ( | const tensor< T > & | t, |
| std::size_t | axis, | ||
| const Order & | larger | ||
| ) |
The largest (larger is std::greater) or smallest elements along axis. Along an inner axis each line is contiguous (extreme); along an outer one the slabs are compared element by element, so the loop runs over contiguous memory rather than down a column.
Definition at line 471 of file tensor_ops.hpp.
|
private |
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).
Each as LAPACK's dlarfg (zlarfg) forms it: H = I − τ w wᴴ, w₀ = 1, with Hᴴ x = β e₁ for a real β of the sign opposite to Re x₀; the identity (τ = 0, an empty reflector) where nothing lies below the diagonal and x₀ is real. Exact inputs stay exact where LAPACK's do (e₃ reflected onto e₁: w = (1, 0, 1), τ = 1; a unit v there holds 1/√2, rounded), and the factors are LAPACK's to rounding. β goes on the diagonal and w₁… below it, where dgeqrf leaves them.
Definition at line 2159 of file decomposition.hpp.
|
private |
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).
Definition at line 2258 of file decomposition.hpp.
|
constexprnoexcept |
The elements one index of parallel_zero or parallel_copy takes of count: zero_chunk, or, from the size at which a buffer takes Linux's 2 MiB pages (16 MiB), one such page's worth, so that no two threads fault in the same page: they waited on each other's fault there (six x86-64 AVX2 cores, 32 MiB zeroed into fresh memory 3.6 → 2.5 ms; the LU and Cholesky of order 2000, their operand copied by eight 4 KiB pages a thread, took 2 to 4 % longer than with one thread's copy).
Definition at line 654 of file parallel.hpp.
|
private |
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.
| [in] | left | The left reflector of each step (acting from its row), empty where none. |
| [in] | right | The right reflector of each step (acting from the column past it), likewise. |
| [in] | left_phase | The unit phase each row took. |
| [in] | right_phase | The unit phase each column past the diagonal took. |
| [in] | n | The order. |
| [out] | ut | U_bᵀ, row-major: its rows are the columns of U_b. |
| [out] | vt | V_bᵀ, likewise. |
Definition at line 4770 of file decomposition.hpp.
|
private |
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.
Definition at line 2363 of file decomposition.hpp.
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).
Definition at line 273 of file decomposition.hpp.
|
private |
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.
Definition at line 1979 of file decomposition.hpp.
| void numerical::detail::gemm | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
C += A·B, every element of C summed in ascending order of the depth.
Products of fewer than gemm_packed_threshold multiply-adds, or shallower than gemm_min_depth, run gemm_rows. The rest run packed (gemm_packed). A double product runs on the AMX where it is in use and suits the product (amx_suits), else on the NEON or AVX tiles: each element a chain of fused multiply-adds; a float one on the NEON tile on arm64's NEON and AMX paths. On the portable path either runs on the vector tile, which rounds c + a·b as the compiler does: each product before its sum with contraction off, as the plain loop then does; fused under Clang's -ffp-contract=on, as the plain loop then is; under GCC's, apart, where the plain loop's scalar statement is fused. Every element is summed in ascending order of the depth on any number of threads. Complex numbers on the x86-64 paths run on their complex tiles (avx2_full_tile_c64, avx512_full_tile_c64 and their _c32 for float); past the portable path on arm64, through the real kernel (gemm_split). Either is the same on any number of threads too, but no longer the plain loop's to the bit: so past the portable path a complex product of any shape runs there, and each element of C is the same whatever the shape of the product it falls in. A caller that splits its columns by the threads (an LU inverse's slices) then has the same bits on any number of them, where a slice narrow enough for the plain loop rounded its complex products otherwise.
| [in] | m,n,k | The shapes: A is m×k, B k×n, C m×n, all row-major. |
| [in] | a,lda,b,ldb | The operands and the distance between their rows. |
| [in,out] | c,ldc | The result and the distance between its rows (C may be a block of a larger matrix that A and B do not overlap). A zero m, n or k leaves C as it is. |
| void numerical::detail::gemm_at | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | at, | ||
| std::size_t | ldat, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
gemm_layouts with A given as Aᵀ (at, k × m, ldat apart).
| void numerical::detail::gemm_atbt | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | at, | ||
| std::size_t | ldat, | ||
| const T * | bt, | ||
| std::size_t | ldbt, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
gemm_layouts with both operands given as their transposes.
| void numerical::detail::gemm_bt | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | bt, | ||
| std::size_t | ldbt, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
gemm_layouts with B given as Bᵀ (bt, n × k, ldbt apart).
| bool numerical::detail::gemm_double_elsewhere | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const double * | a, | ||
| std::size_t | lda, | ||
| const double * | b, | ||
| std::size_t | ldb, | ||
| double * | c, | ||
| std::size_t | ldc | ||
| ) |
A double product on the path in use (numerical::gemm_path_in_use), unless it is the portable one, B and A laid out as B and A say (b_rows, a_rows and their transposes).
|
inline |
gemm_double_elsewhere of a row-major B.
|
inline |
gemm_float_elsewhere of row-major operands.
| void numerical::detail::gemm_kernel | ( | std::size_t | kc, |
| const T * | ap, | ||
| const T * | bp, | ||
| std::size_t | ldb, | ||
| gemm_tile< T > & | tile | ||
| ) |
| void numerical::detail::gemm_large | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
A product past the plain loop: on the AMX or x86-64 tiles where gemm_elsewhere takes it, else on the vector tile.
gemm_large_complex of row-major operands.
| void numerical::detail::gemm_large_layouts | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const double * | a, | ||
| std::size_t | lda, | ||
| const double * | b, | ||
| std::size_t | ldb, | ||
| double * | c, | ||
| std::size_t | ldc | ||
| ) |
A double product past the plain loop, B and A laid out as B and A say (gemm_layouts).
| void numerical::detail::gemm_large_layouts | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const float * | a, | ||
| std::size_t | lda, | ||
| const float * | b, | ||
| std::size_t | ldb, | ||
| float * | c, | ||
| std::size_t | ldc | ||
| ) |
A float product past the plain loop, B and A laid out as B and A say (gemm_layouts).
| void numerical::detail::gemm_large_layouts | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
A complex product past the plain loop, B and A laid out as B and A say (gemm_layouts).
| void numerical::detail::gemm_layouts | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
gemm with either operand given as its transpose (A: a_rows or a_transposed, B: b_rows or b_transposed), so that a product by a transpose (a syrk's P Pᵀ, BLAS's transa and transb, a NumPy array in Fortran order) needs no copy of it. The tiles are gemm's, the operands packed (or, through gemm_split, copied into their parts) from where they lie, and so are the bits; a product gemm runs in the plain loop copies a transposed operand row-major and runs there, small as it is.
| void numerical::detail::gemm_pack_at | ( | const T * | at, |
| std::size_t | ldat, | ||
| std::size_t | rows, | ||
| std::size_t | kc, | ||
| T * | out | ||
| ) |
The slivers of Mr rows of A (rows × kc), column by column, from its transpose at (kc × rows row-major, ldat apart): each column of a sliver is a run of a row of at; a last partial sliver is padded with zeros, as gemm_pack_a lays them out.
| void numerical::detail::gemm_pack_bt | ( | const T * | bt, |
| std::size_t | ldbt, | ||
| std::size_t | kc, | ||
| std::size_t | cols, | ||
| T * | out | ||
| ) |
gemm_pack_b of B = Bᵀᵀ, from the row-major bt (cols × kc, leading dimension ldbt): each column of a sliver is a row of bt, read along it; a last partial sliver is padded with zeros.
| void numerical::detail::gemm_packed | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc, | ||
| const T * | a_packed = nullptr |
||
| ) |
The packed product with the tile Kernel computes: for each block of B's columns and each panel of the depth, both operands are packed (in parallel, by the blocks the product takes), then tasks of gemm_mc rows by gemm_slivers_per_task slivers of columns run in parallel, each writing a disjoint block of C, one tile at a time. Given a_packed (A as pack_left packs it for Kernel, a depth of one panel), A is read from there and only B is packed: the same tiles, the same bits. B is B's layout (b_rows, b_transposed), A A's (a_rows, a_transposed): the packing alone differs, so the bits do not.
| void numerical::detail::gemm_packed_triangle | ( | std::size_t | n, |
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc, | ||
| bool | lower | ||
| ) |
gemm_packed restricted to a triangle of the square C = A·B (n × n): its lower one (j ≤ i) where lower, else its upper one (j ≥ i); the other is neither written nor depended on. The same packing (A's row blocks with no tile in the triangle left out) and the same tiles, so each element's sum is the product's; a tile across the diagonal runs on a copy of its block of C, the triangle's part copied back. A syrk takes the triangle in one pass, its operands packed once a panel, where products by halves packed them again for every block.
|
inlinenoexcept |
The most capable path products may take, as numerical::limit_gemm_path set it: none is probed until a product (or a question) reaches it, so a cap below the AMX keeps its instructions from ever running.
|
inlinenoexcept |
The path the next double product takes: the most capable one under the cap that runs here (each path's probe is made once). Kept with its cap, one load, where each call walked the paths down from the cap through their probes' guards: every small step of a decomposition asks (an LU of order 100 through dgetrf_ on one AVX2 core spent 5 % of its time here). A race stores the same value twice.
|
inlinenoexcept |
The path gemm_path_current found last, packed with the cap it was found under (cap · 8 + path); 0xFF before any.
| void numerical::detail::gemm_rows | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
C += A·B by the plain i-k-j loop (gemm_rows_rounded), each element rounded as the tile of the path in use rounds it, so that it is the same whatever the shape of the product it falls in: a caller that splits its columns by the threads (an LU inverse) then has the same bits on any number of them.
| void numerical::detail::gemm_rows_rounded | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| const T * | b, | ||
| std::size_t | ldb, | ||
| T * | c, | ||
| std::size_t | ldc | ||
| ) |
C += A·B by the plain i-k-j loop, parallel over the rows of C. Past a depth of 1 (narrow_least_depth), the rows run narrow_group at a time by blocks (narrow_rows), the last rows one by one: the same sums, six to ten times as fast from order 4 to 16 (M1 Pro) as the loop of unknown width, whose sums the compiler cannot keep in registers. A depth of 1 streams as fast through that loop.
| [in] | m,n,k | The shapes: A is m×k, B k×n, C m×n, all row-major. |
| [in] | a,lda,b,ldb | The operands and the distance between their rows. |
| [in,out] | c,ldc | The result and the distance between its rows. Each term rounded as Mode says (term_rounding). |
| void numerical::detail::gemm_split | ( | std::size_t | m, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| const Z * | a, | ||
| std::size_t | lda, | ||
| const Z * | b, | ||
| std::size_t | ldb, | ||
| Z * | c, | ||
| std::size_t | ldc | ||
| ) |
C += A·B for complex numbers through the real kernel, by blocks of split_columns columns and split_depth of depth, on the real path's units and threads (a scalar complex tile has none). Within a block, each element of C takes the products of the real parts of A in ascending order of the depth, then those of the imaginary parts: by four real products of the block's parts (split_parts_block) or one of its expanded operands (split_expanded_block), as split_by_parts chooses, the same bits either way. A and B laid out as A and B say, read where they lie as the parts are copied.
| void numerical::detail::gemm_triangle | ( | std::size_t | n, |
| std::size_t | k, | ||
| const double * | a, | ||
| std::size_t | lda, | ||
| const double * | b, | ||
| std::size_t | ldb, | ||
| double * | c, | ||
| std::size_t | ldc, | ||
| bool | lower | ||
| ) |
The triangle of the double product C = A·B (gemm_packed_triangle) on the tile of the path in use, A and B laid out as A and B say.
|
inline |
|
inline |
|
inlinenoexcept |
The threads the pool runs, the caller's included (at least 1), known without starting the pool: one per physical core where the system tells them (else one per hardware thread). Two threads on one core's hyperthreads share its multiply-add units: a product of order 128 on eight AVX2 cores with sixteen hyperthreads took 140 µs on sixteen threads, 71 on eight; order 2048, the same either way. Asked once, as the pool's workers are started once (glibc reads a file for the hardware threads on every call).
Definition at line 210 of file parallel.hpp.
The real symmetric or complex Hermitian matrix whose lower triangle is the square a's, and of its diagonal the real part.
Definition at line 6653 of file decomposition.hpp.
| 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).
| [in] | block | The block's first element; its rows n apart. |
| [in] | n | The row stride. |
| [in] | m | The block's order. |
| [in] | v | The vector (m elements). |
| [out] | w | B v (m elements). |
| [in,out] | scratch | The partials, kept by the caller from one step of a reduction to the next: a fresh buffer of their size each step is mapped anew, page by page. |
Definition at line 2879 of file decomposition.hpp.
|
inline |
hermitian_pseudo_inverse in double: V Λ⁺ Vᵀ from the eigendecomposition, cheaper.
Definition at line 6759 of file decomposition.hpp.
|
inline |
hermitian_pseudo_inverse in std::complex<double>: V Λ⁺ Vᴴ from the eigendecomposition, cheaper.
Definition at line 6768 of file decomposition.hpp.
| 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.
Definition at line 6751 of file decomposition.hpp.
|
inline |
The same for a real row, in one pass through kernel.
Definition at line 2835 of file decomposition.hpp.
|
inline |
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.
Definition at line 2803 of file decomposition.hpp.
| 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.
Definition at line 2740 of file decomposition.hpp.
hermitian_singular_values in double: the magnitudes of the eigenvalues, cheaper.
Definition at line 6735 of file decomposition.hpp.
|
inline |
hermitian_singular_values in std::complex<double>: the magnitudes of the eigenvalues, cheaper.
Definition at line 6741 of file decomposition.hpp.
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).
Definition at line 6716 of file decomposition.hpp.
| buffer< T > numerical::detail::hessenberg | ( | buffer< T > & | a, |
| std::size_t | n, | ||
| std::size_t | low, | ||
| std::size_t | high, | ||
| bool | form | ||
| ) |
hessenberg_blocked from an active order of hessenberg_crossover, else hessenberg_unblocked.
| buffer< T > numerical::detail::hessenberg_blocked | ( | buffer< T > & | a, |
| std::size_t | n, | ||
| std::size_t | low, | ||
| std::size_t | high, | ||
| bool | form | ||
| ) |
hessenberg_unblocked by panels of hessenberg_block columns, as LAPACK's dgehrd and dlahr2: a panel's reflectors are formed one column at a time, each column first updated by those before it through Y = A V T; the rest of the matrix then takes the panel from both sides through the matrix product kernel, A ← (I − V Tᴴ Vᴴ)(A − Y Vᴴ). Its arguments and answer are hessenberg_unblocked's.
| buffer< T > numerical::detail::hessenberg_unblocked | ( | buffer< T > & | a, |
| std::size_t | n, | ||
| std::size_t | low, | ||
| std::size_t | high, | ||
| bool | form | ||
| ) |
Reduces the active block low … high − 1 of the balanced row-major n × n a to upper Hessenberg form by Householder reflections one column at a time, as LAPACK's dgehd2: A = Q H Qᴴ, the entries below the subdiagonal set to zero.
| [in,out] | a | The matrix; on return H. |
| [in] | n | Its order. |
| [in] | low | The active block's first row and column. |
| [in] | high | One past its last. |
| [in] | form | Whether to form Q. |
Qᵀ, row-major n × n (the identity outside the active block), or empty without form: a QR sweep updates columns of Q, which are then contiguous. The unit Householder vector v with (I − 2 v vᵀ) x = β e₁, in place of x, for a real x.
β = −sign(x₀) ‖x‖. Definition at line 4841 of file decomposition.hpp.
|
private |
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.
Blocked as LAPACK's dgeqrf: the reflectors of a panel of qr_block columns act on the panel one by one, then on the columns right of it at once, as I − V T Vᴴ through the matrix product kernel; Q is formed block by block the same way.
| [in] | a | The matrix, row-major. |
| [in] | rows | Its row count. |
| [in] | cols | Its column count. |
| [in] | complete | Whether Q is the whole rows × rows product of the reflections, its columns past k an orthonormal basis of the complement of A's column space (numpy.linalg.qr(mode="complete")). |
Q and R. Definition at line 2419 of file decomposition.hpp.
|
private |
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).
Definition at line 2329 of file decomposition.hpp.
| 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.
Definition at line 6675 of file decomposition.hpp.
Whether x is finite: both of its parts.
| [in] | x | The entry. |
Definition at line 473 of file decomposition.hpp.
Is ch an ASCII letter (a valid index label)?
Definition at line 341 of file tensor.hpp.
A NaN? Never for a scalar without one.
Definition at line 331 of file tensor_ops.hpp.
A NaN?
Definition at line 337 of file tensor_ops.hpp.
The Frobenius norm of the strict off-diagonal of the n × n a.
Definition at line 4262 of file decomposition.hpp.
|
inline |
One Jacobi rotation of the n × n a zeroing a[p][q], accumulated into v.
Definition at line 4275 of file decomposition.hpp.
|
inlineprivate |
Joins the halves of nodes' block index in work, then lets them go.
Definition at line 5873 of file decomposition.hpp.
|
inlineprivate |
Joins the two decomposed halves of node, one and two, by merge_arrow, node started (start_bidiagonal_node).
Definition at line 5794 of file decomposition.hpp.
|
inline |
Joins the two solved halves of node by merge_rank_one, in work.
Definition at line 4082 of file decomposition.hpp.
|
noexcept |
|
inlinenoexcept |
left_kernel_for a double product: past gemm_rows, as gemm_elsewhere chooses, of a depth of one panel; plain otherwise.
|
noexcept |
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.
Definition at line 598 of file decomposition.hpp.
|
private |
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.
| [in] | a | The m × n matrix A. |
| [in] | b | The right-hand sides, one per column. |
| [in] | cutoff | Where a singular value counts as zero. |
| [out] | r | Where given, the square R factorised, of A divided by 2^exponent, and that exponent: its singular values are A's over 2^exponent. |
Definition at line 6334 of file decomposition.hpp.
|
private |
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.
Definition at line 6426 of file decomposition.hpp.
|
constexprnoexcept |
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.
Definition at line 127 of file decomposition.hpp.
| void numerical::detail::matrix_times_vector | ( | const T * | a, |
| std::size_t | rows, | ||
| std::size_t | cols, | ||
| const T * | x, | ||
| T * | y | ||
| ) |
y = A x for the row-major rows × cols a: each element the dot product of its row with x, as chunked_dot sums it (a row of a chunk or less is one dot_product), the rows in parallel by groups (row_dots, which share each load of x). The same on any number of threads. A double matrix of amx_dot_least elements to fewer than amx_dot_most, on the AMX path, runs on the AMX (amx_dot_rows_f64): each element then its eight lanes of fused multiply-adds, added in order.
|
inlineprivate |
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.
Deflation first (dlasd2): a negligible zⱼ leaves dⱼ and its vectors as they are, a dⱼ near zero is rotated into the zero pole, and two nearly equal dⱼ are rotated together so that one of them deflates. The remaining roots are found by singular_secular_root, in parallel; z is recomputed from them (Löwner's formula), so that the vectors are orthogonal to working precision, and they reach u and v by a matrix product each.
| [in,out] | d | The diagonal; on return, the singular values. |
| [in] | z | The first row. |
| [in] | alpha,beta | The entries the arrow was built from: they set the deflation tolerance. |
| [in,out] | u | The left basis; on return, the left singular vectors. |
| [in,out] | v | The right basis; on return, the right singular vectors. |
| [in] | n | The arrow's order. |
| [in] | m | The right basis's order (n, or n + 1 with a null column). |
| [in] | split | The row of u the arrow's first row is (its upper half above, its lower below): the bases are block diagonal by it, and each half's rows take only the columns that reach them (dlasd3 multiplies the halves apart the same way). |
| [in,out] | work | The buffers it works in, the caller's thread's (thread_arrow_workspace). |
Definition at line 5448 of file decomposition.hpp.
|
inline |
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.
Deflation first (dlaed2): a negligible zⱼ leaves dⱼ and its column as they are, and two nearly equal poles are rotated so that one of them does. The remaining k roots are found by secular_root, in parallel; z is then recomputed from them (Löwner's formula, Gu and Eisenstat), so that the eigenvectors (dⱼ − λᵢ)⁻¹ ẑⱼ are orthogonal to working precision, and they reach q by a matrix product per half: q is block diagonal, its first split columns zero below row split and the others above it, and only a deflating rotation mixes a column of each (dlaed3 multiplies the halves apart the same way).
Definition at line 3790 of file decomposition.hpp.
|
inline |
zt ← left · zt, n × n, row-major, for a real left factor.
Definition at line 4205 of file decomposition.hpp.
|
inline |
zt ← left · zt for a complex zt: the real and imaginary parts by a real product each.
Definition at line 4215 of file decomposition.hpp.
| 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.
Definition at line 4191 of file decomposition.hpp.
|
constexpr |
|
constexpr |
The order from which a real Schur form is taken by multishift_schur, below which real_schur is faster (eight arm64 cores: order 128 1.85 against 2.60 ms; 224, 8.9 against 8.6; with the vectors, 192, 8.1 against 7.1).
| void numerical::detail::multishift_schur | ( | const schur_target< T > & | target, |
| std::complex< double > * | w, | ||
| bool | recursive, | ||
| lapack_index | steps | ||
| ) |
The Schur form of target's active block by the multishift QR algorithm with aggressive early deflation, as LAPACK's dlaqr0 and zlaqr0 (dlaqr4 without recursive): each iteration deflates what it can from a window at the bottom of the active block (early_deflation) and, unless that deflated enough, sweeps the block with the window's other eigenvalues as shifts (chase_bulges). The eigenvalues are left in w (one per row; a real matrix's complex pairs adjacent, the positive imaginary part first).
| [in] | target | The Hessenberg matrix, its active block, and the Schur vectors to update. |
| [out] | w | The eigenvalues, one per row of the matrix. |
| [in] | recursive | Whether a wide deflation window is itself taken by these iterations (dlaqr0), or by the single-shift ones (dlaqr4). |
| [in] | steps | The iterations allowed, LAPACK's 30 · max(10, order) where zero. |
| convergence_error | If the iterations did not converge within steps. |
|
noexcept |
C += A·B for R rows and N columns of C, each sum in ascending order of the depth, as the plain loop's: the block held in registers, of a size the compiler knows; each term rounded as Mode says. Out of line: its unrolled sizes, inlined, spent GCC's growth budget, and the LU factorisation, which runs none of them, lost 3 to 8 % at every order (one arm64 thread, GCC 14).
|
noexcept |
C += A·B for R whole rows of C of n columns: blocks of narrow_columns, the last narrower.
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.
Definition at line 350 of file decomposition.hpp.
|
inlineprivate |
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).
Definition at line 6252 of file decomposition.hpp.
| C numerical::detail::numpy_divide | ( | const C & | a, |
| const C & | b, | ||
| bool | fused | ||
| ) |
a / b of complex numbers as NumPy's array loop computes it: Smith's division by the larger part of b (an exact zero divided part by part, to an infinity or a NaN), each x + y·z of it one fused multiply-add where fused (NumPy's build contracts them on arm64, not on x86-64).
Definition at line 132 of file tensor_ops.hpp.
|
inlinenoexcept |
Whether NumPy's complex product fuses its multiply-adds on this processor: on arm64, and on x86-64 where its AVX2 loop runs (FMA with it); elsewhere its loop rounds each product.
Definition at line 154 of file tensor_ops.hpp.
|
inlinenoexcept |
Whether NumPy's complex quotient fuses its multiply-adds, as its PyPI wheels are built: on macOS arm64, whose clang contracts them (FMA is arm64's baseline); not on Linux, whose GCC does not, nor on x86-64, whose baseline has no FMA.
Definition at line 180 of file tensor_ops.hpp.
|
inlinenoexcept |
Whether NumPy's complex product fuses its multiply-adds in its scalar loop (where numpy_loop_is_scalar holds), as its PyPI wheels are built: that loop sits beside the vector one, and clang (macOS) and MSVC (Windows) contract it where the vector loop has FMA; GCC (Linux, ISO C) does not.
Definition at line 168 of file tensor_ops.hpp.
|
noexcept |
Whether NumPy multiplies a by b in its scalar loop rather than its vector one: when the result is one element and an operand of rank above 0 is broadcast to a higher rank, NumPy's iterator hands its loop an output step of 0, which the vector loop refuses.
Definition at line 195 of file tensor_ops.hpp.
x · y + z: one fused multiply-add where fused, else the product rounded first (statements apart, so no compiler contracts them).
Definition at line 96 of file tensor_ops.hpp.
| C numerical::detail::numpy_multiply | ( | const C & | a, |
| const C & | b, | ||
| bool | fused | ||
| ) |
a · b of complex numbers as NumPy's array loop computes it: ai·bi and ai·br rounded, then ar·br − ai·bi and ar·bi + ai·br, each by a fused multiply-add where fused (NumPy's vectors on arm64 and on x86-64 with FMA), else each product rounded first. No recovery of an infinity from a NaN, as std::complex's (C Annex G) makes: 1j · (inf + inf·j) is NaN, (1e308 + 1e308j) · (1e308 − 1e308j) inf + inf·j when fused.
Definition at line 116 of file tensor_ops.hpp.
|
private |
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.
| [in,out] | w | The matrix whose columns are orthogonalised. |
| [in,out] | v | The rotations, accumulated (start from the identity). |
| [in] | n | The order. |
| [in] | max_sweeps | The bound of sweeps. |
| convergence_error | If a sweep still rotated at the bound. |
Definition at line 4648 of file decomposition.hpp.
|
inline |
The count entries stride apart at v as a real_operand.
|
inline |
The same as a complex_operand.
| bool numerical::detail::operator== | ( | const std::vector< T, default_init_allocator< T > > & | buffer, |
| const std::vector< T > & | other | ||
| ) |
Element-wise equality of a numerical::buffer and a std::vector, found by argument-dependent lookup through the allocator: code that compared a tensor's elements with a std::vector keeps doing so.
Definition at line 144 of file buffer.hpp.
a · b, as std::operator* computes it, in line. libc++'s complex product (the four products and their sums as statements of their own, so unfused) is copied here, its rescue of a product whose parts are both NaN left to it: where it was a call, out of line in a translation unit as large as the Python extension, it took 45 % of a complex LU of order 48. Elsewhere (libstdc++, whose product the compiler expands) and for a real T, the operator itself.
| 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.
Definition at line 82 of file decomposition.hpp.
x over the pivot whose pivot_divisor is divisor: the same quotient wherever the LU takes it.
Definition at line 74 of file decomposition.hpp.
| packed_left< T > numerical::detail::pack_left | ( | std::size_t | m, |
| std::size_t | , | ||
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| buffer< T > & | |||
| ) |
The m × k a packed once for products with it of about n columns each (gemm_with), in storage (grown as needed, kept by the caller from one call to the next): on the tile gemm would choose for m × n × k. A scalar other than double (here), a depth of more than one panel, or a product gemm_rows takes, leaves A to be read in place.
|
inline |
pack_left for double, on the tile left_kernel_for chooses.
| const T * numerical::detail::pack_left_for | ( | std::size_t | m, |
| std::size_t | k, | ||
| const T * | a, | ||
| std::size_t | lda, | ||
| buffer< T > & | storage | ||
| ) |
A packed for Kernel into storage, in parallel by its blocks of rows, as gemm_packed packs it.
|
inline |
|
inline |
|
inlineprivate |
The labels of labels[i] and labels[j] that the output or another operand carries, in their order of appearance, once each: what their contraction keeps.
Definition at line 700 of file tensor.hpp.
|
inline |
|
inline |
|
private |
Sum of count values from first at stride stride, by NumPy's pairwise summation: below 8 values in order; up to 128 in eight interleaved partial sums, added pairwise, then the remainder in order; above, the two halves (the first a multiple of 8 long) summed apart.
Definition at line 216 of file tensor.hpp.
|
private |
pairwise_sum of width neighbouring columns at once, into out: column q's sum that of pairwise_sum(first + q, count, stride), to the bit, each step adding a whole row of the columns (contiguous, stride apart) where a column was walked down stride at a time.
Definition at line 257 of file tensor.hpp.
|
private |
The pairwise sums of outer blocks of count rows of inner elements (row-major), each column of a block summed down its rows, into out (outer × inner): one column at a time where inner is 1, else by rows (pairwise_sum_rows), the same bits.
Definition at line 322 of file tensor.hpp.
Copies the count elements at in to out, split across the threads the work pays for, as parallel_zero does: a fresh copy's pages are faulted in by several threads, where a sequential copy pays every fault on one (a 2 MiB copy into fresh memory took 0.85 ms on one x86-64 core in a container, a third of the Cholesky factorisation of order 512).
Definition at line 703 of file parallel.hpp.
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).
Definition at line 545 of file decomposition.hpp.
| void numerical::detail::parallel_for | ( | std::size_t | count, |
| std::size_t | work_per_index, | ||
| const Body & | body | ||
| ) |
Runs body(i) for every i in [0, count).
The range is partitioned across as many threads as its work pays for (threads_for), and runs sequentially when that is fewer than two. body must be safe to invoke concurrently for distinct i (it is when each call touches a disjoint part of the output).
An exception body throws in a worker is kept and rethrown here once every worker has joined (the first chunk's, in index order): one escaping a worker thread would call std::terminate. A worker the system cannot start is a std::system_error, thrown once the others have joined.
| [in] | count | The number of iterations. |
| [in] | work_per_index | The work of one iteration, in multiply-adds. |
| [in] | body | The per-index work (void(std::size_t)). |
Definition at line 580 of file parallel.hpp.
| void numerical::detail::parallel_for_beside | ( | std::size_t | count, |
| std::size_t | work_per_index, | ||
| const Body & | body, | ||
| const Side & | side | ||
| ) |
Runs side on the caller while the pool's workers run body(i) for every i in [0, count), then lets the caller take the chunks left: one task beside a parallel loop, each touching what the other does not (a factorisation's next panel beside its trailing update). side's own parallel loops find the pool busy and run in its thread. Where the work pays for fewer than two threads, or the pool is busy, side runs, then the loop, on the caller. An exception of side comes back after the loop's chunks have joined, before theirs.
| [in] | count | The number of iterations. |
| [in] | work_per_index | The work of one iteration, in multiply-adds. |
| [in] | body | The per-index work (void(std::size_t)). |
| [in] | side | The caller's task (void()). |
Definition at line 618 of file parallel.hpp.
Sets the count elements at out to zero, split across the threads the work pays for: a fresh result's pages are first touched, and faulted in, by several threads, where a sequential fill pays every fault on one (a product on eight arm64 threads: 16 % faster at order 512, 7 % at 1024). The zero is a constant in each chunk's fill, which the compiler turns into the C library's zeroing (on arm64 it clears whole cache lines without reading them; a value passed in, stored element by element, was 4 % slower at order 1024 on one thread).
Definition at line 682 of file parallel.hpp.
|
private |
The elements of data, whose axes carry the labels from (extents in sizes), laid out row-major along the labels to, a permutation of from that differs from it (so at least one label): an odometer over to's outer labels, a strided sweep of its last.
Definition at line 541 of file tensor.hpp.
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.
Definition at line 320 of file decomposition.hpp.
|
inlinenoexcept |
The physical cores this process may run on, or 0 where the system does not say. Linux: core_sets; macOS: hw.physicalcpu.
Definition at line 180 of file parallel.hpp.
pivot_divisor of a complex pivot: its reciprocal.
Definition at line 67 of file decomposition.hpp.
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.
Definition at line 60 of file decomposition.hpp.
|
inlinenoexcept |
Whether a double product of m rows, n columns and depth k, past gemm_rows's threshold, still runs faster in the plain loop than packed on the path in use: packing B's slivers for so few rows of A, or so short a depth, did not pay. One core, the plain loop's time over the packed product's: the NEON tile at a depth of 16 or less whatever the shape (0.5 to 0.6 up to 2000 × 2000), at 64 or less for 64 rows or fewer (0.4 to 0.9; the same under GCC 14 and on eight threads); the AVX2 tile for 16 rows or fewer to a depth of 64, or 32 rows to 32 (0.2 to 0.8, six cores too), as wide as plain_loop_least_width. The AMX's packing pays from 16 rows.
|
noexcept |
plain_loop_pays for a product of T: double's alone, the one it was measured for.
|
inlinenoexcept |
plain_loop_pays for double.
√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.
Definition at line 378 of file decomposition.hpp.
|
inline |
The Schur form of t's active block by real_schur, its eigenvalues to w (one per row of t).
|
inline |
The same by complex_schur.
|
inlineprivate |
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.
nodes. Definition at line 5687 of file decomposition.hpp.
|
inline |
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.
Definition at line 4012 of file decomposition.hpp.
|
constexprnoexcept |
The pool's threads from the logical hardware threads and the physical cores (0: unknown): the cores where known, at least 1.
Definition at line 196 of file parallel.hpp.
Nothing, for a tile that prefetches no C (prefetches_next_row_of_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).
| [in] | a | A. |
| [in] | b | B. |
| [in] | c | C. |
η. Definition at line 3586 of file decomposition.hpp.
|
inlineprivate |
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.
Definition at line 2459 of file decomposition.hpp.
|
noexcept |
The real part of a complex scalar.
Definition at line 287 of file decomposition.hpp.
The real part of a real (or general) scalar as a double.
Definition at line 280 of file decomposition.hpp.
|
inline |
The real Schur form of the upper Hessenberg t's active block, as LAPACK's dlahqr: Francis's implicit double-shift QR, with the deflation criterion of Ahues and Tisseur and exceptional shifts after 10 and 20 iterations without deflation. The eigenvalues are left in re and im, a complex pair's in two consecutive places, the positive imaginary part first.
| [in] | t | The Hessenberg matrix, its active block, and the Schur vectors to update. |
| [in,out] | re | The eigenvalues' real parts (the isolated ones already in place). |
| [out] | im | Their imaginary parts. |
| [in] | steps | The iterations allowed between two deflations. |
| convergence_error | If an eigenvalue has not converged within steps. |
A real subdiagonal entry of the tridiagonal: kept with its sign (the rotations take any sign).
| [in] | t | The entry. |
t. Definition at line 3129 of file decomposition.hpp.
|
inline |
A complex subdiagonal entry made real by the next row's unit phase: next = delta · conj(phase(t)).
| [in] | t | The entry. |
| [in] | delta | The phase of the row above. |
| [out] | next | The phase of the row below. |
|t|. Definition at line 3143 of file decomposition.hpp.
|
inline |
The eigenvectors of the real quasi-triangular t (row-major n × n, its 2 × 2 blocks standard) as the rows of Xᵀ, as LAPACK's dtrevc: real for a real eigenvalue, solved in real arithmetic; for a pair re ± i im (im [k] > 0 the first) the first's vector, the second's its conjugate.
|
inline |
|
inline |
| tensor< T > numerical::detail::reduce_along | ( | const tensor< T > & | t, |
| std::size_t | axis, | ||
| const Reduce & | reduce | ||
| ) |
Reduces each line of t along axis by reduce, which takes the line's first element, its length and its stride.
Definition at line 300 of file tensor_ops.hpp.
|
private |
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.
Definition at line 2127 of file decomposition.hpp.
|
private |
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).
Definition at line 2209 of file decomposition.hpp.
|
private |
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).
Definition at line 2292 of file decomposition.hpp.
|
private |
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.
Row by row: the projection w = vᴴ A accumulates one contiguous row at a time, then each row takes −2 vᵢ w. A column-by-column walk strides by cols on every access and runs several times slower. The columns are independent, so chunks of them run in parallel, each element's arithmetic in the same order.
| [in,out] | a | The array, cols entries per row. |
| [in] | cols | Its row length. |
| [in] | v | The unit reflector. |
| [in] | first | The first row it acts on. |
| [in] | from | The first column it acts on. |
Definition at line 1937 of file decomposition.hpp.
| void numerical::detail::reflect_three_rows | ( | T * | r1, |
| T * | r2, | ||
| T * | r3, | ||
| std::size_t | from, | ||
| std::size_t | to, | ||
| const T & | tau, | ||
| const T & | v2, | ||
| const T & | v3, | ||
| bool | adjoint | ||
| ) |
A reflection I − τ v vᴴ, v = (1, v₂, v₃), on the entries from … to − 1 of three rows: its adjoint from the left with adjoint (rows of H), else from the right (rows of Zᵀ, which are columns of Z). A reflection of two passes r3 = r2 and v3 = 0. One loop over the entries, which vectorises.
| void numerical::detail::reflect_two_rows | ( | T * | r1, |
| T * | r2, | ||
| std::size_t | from, | ||
| std::size_t | to, | ||
| const T & | tau, | ||
| const T & | v2, | ||
| bool | adjoint | ||
| ) |
reflect_three_rows for a reflection of two entries, v = (1, v₂), on two rows: its own loop, which vectorises where the third row aliasing the second would not.
|
inlineprivate |
Releases the buffers of work longer than kept_workspace.
Definition at line 5409 of file decomposition.hpp.
|
inline |
Releases the buffers of work longer than kept_workspace.
Definition at line 3768 of file decomposition.hpp.
|
inline |
Moves the diagonal block of the real quasi-triangular t at ifst to ilst by exchanges of adjacent blocks, the rows of qt (Qᵀ) turning with it, as LAPACK's dtrexc. ilst becomes where the block ended: short of its goal where an exchange was refused (then.
|
inline |
require_finite for double, its elements contiguous: all_finite.
Definition at line 577 of file decomposition.hpp.
|
inline |
The same for std::complex<double>, whose array holds its parts in order ([complex.numbers]).
Definition at line 586 of file decomposition.hpp.
The one non-finite policy, checked where every decomposition starts.
| [in] | a | The input. |
| [in] | what | The decomposition, for the message. |
| non_finite_error | If an entry is infinite or NaN. |
Definition at line 483 of file decomposition.hpp.
|
inline |
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.
Definition at line 631 of file decomposition.hpp.
|
inline |
The same for std::complex<double>, whose array holds its parts in order ([complex.numbers]).
Definition at line 640 of file decomposition.hpp.
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).
| [in] | a | The square input. |
| [in] | what | The decomposition, for the message. |
| non_finite_error | If an entry read is infinite or NaN. |
Definition at line 618 of file decomposition.hpp.
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.
float). Definition at line 334 of file decomposition.hpp.
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.
Definition at line 308 of file decomposition.hpp.
|
inlineprivate |
Columns a and b of the row-major m (cols per row) rotated: a ← c a + s b, b ← c b − s a.
Definition at line 5366 of file decomposition.hpp.
√(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]].
Definition at line 426 of file decomposition.hpp.
|
inline |
The operations split into runs of chained rotations, every other operation a run of one.
Definition at line 3373 of file decomposition.hpp.
|
inlinenoexcept |
rounded_apart_add of complex numbers: the product in line (operator_times, the operator's own rounding, unfused), where it was a call that took a sixth of a complex LU.
sum + a·b as the vector tile rounds it (vector_tile_kernel), so that the plain loop's elements are its own: under Clang's -ffp-contract=on one expression, fused in both; under GCC's, whose vector extensions fuse nothing while a scalar statement's multiply-add fuses, the product apart, a statement of its own.
|
inlinenoexcept |
sum + a·b as the compiler rounds it (rounded_apart_add).
|
inlinenoexcept |
|
inline |
row_dots for double: on the x86-64 units each row by their dot_product (wider than four rows of the portable vectors), elsewhere four rows at a time.
|
inline |
row_dots for std::complex<double>, through complex_row_dots.
|
inline |
row_dots for std::complex<float>, through complex_row_dots.
| void numerical::detail::row_dots | ( | const T * | m, |
| std::size_t | ld, | ||
| std::size_t | rows, | ||
| const T * | x, | ||
| std::size_t | count, | ||
| T * | out | ||
| ) |
dot_product of each row of m with x, for a type no native vector holds.
|
inlinenoexcept |
Whether the tile a double or float product past the plain loop runs on fuses each term (gemm_rows rounds its own so): double on arm64's NEON and AMX paths, float on the NEON tile there. The x86-64 tiles' double takes term_rounding::x86_fused; the vector tile, and so the portable path, rounds as the compiler does, as the plain loop then does. On arm64 every cap but the portable one runs the NEON or the AMX path, so the cap alone answers: one load, where the path in use took an order-4 product on one arm64 core from 10.6 to 13.3 ns (the cap, 11.5).
|
private |
scaled_by_two for a complex number, part by part.
Definition at line 6284 of file decomposition.hpp.
x times two to the exponent, exactly where the result is a normal number (no power of two is formed, so none overflows).
Definition at line 6276 of file decomposition.hpp.
scaled_norm of a whole vector.
Definition at line 414 of file decomposition.hpp.
The same for double: the plain sum of squares in one pass of dot_product where it holds.
Definition at line 397 of file decomposition.hpp.
|
inline |
The same for std::complex<double>: dot_product over the parts, which an array of them holds in order ([complex.numbers]).
Definition at line 405 of file decomposition.hpp.
The Euclidean norm of the count elements at x, by norm_at_scale for this type.
Definition at line 390 of file decomposition.hpp.
|
inline |
The real Schur form of t's active block: real_schur below multishift_from, else multishift_schur, as LAPACK's dhseqr chooses. steps bounds real_schur's iterations.
|
inline |
The complex Schur form likewise: complex_schur or multishift_schur.
|
inline |
Root i of 1/ρ + Σⱼ zⱼ²/(dⱼ − λ) = 0 (d strictly ascending, ‖z‖ = 1), the i-th eigenvalue of diag(d) + ρ z zᵀ, by solve_secular.
| [in] | d | The poles. |
| [in] | z | The weights. |
| [in] | rho | ρ. |
| [in] | i | The root. |
| [in] | max_iterations | The bound on the steps. |
| [out] | delta | dⱼ − λ for every j. |
λ. Definition at line 3695 of file decomposition.hpp.
| T numerical::detail::settled_zero | ( | T | value, |
| const T * | first, | ||
| std::size_t | count, | ||
| std::size_t | stride, | ||
| const Order & | larger | ||
| ) |
The extreme value of count elements stride apart from first, its zero's sign settled as IEEE 754-2019's maximum and minimum settle it, and NumPy's on arm64: no zero is larger than the other, so the first met stayed; max gives +0 where one is among them, min −0, whatever their order. A NaN, or any other value, as it is.
settled_zero of a type with one zero: value.
Definition at line 255 of file tensor_ops.hpp.
|
inline |
A shape as NumPy prints it, (2, 3), for a message.
Definition at line 37 of file tensor_ops.hpp.
|
inlineprivate |
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.
Solved by solve_secular in μ = σ² − d_o², the poles' offsets (dⱼ − d_o)(dⱼ + d_o) computed from their two accurate factors, so that both dⱼ − σ and dⱼ + σ come out with small relative errors.
| [in] | d | The diagonal (the poles). |
| [in] | z | The first row (the weights). |
| [in] | i | The singular value. |
| [in] | max_iterations | The bound on the steps. |
| [out] | below | dⱼ − σ for every j. |
| [out] | above | dⱼ + σ for every j. |
σ. Definition at line 5345 of file decomposition.hpp.
|
constexprnoexcept |
Whether the worker pinned to core worker sits out a run whose caller is on core caller: one thread per core, as two on one core's hyperthreads share its multiply-add units. A free worker never does.
Definition at line 78 of file parallel.hpp.
|
inline |
The solution X (n1 × n2, row-major in x) of TL X − X TR = B, n1 and n2 at most 2, by Gaussian elimination with complete pivoting on its Kronecker form, a pivot below ε times the largest coefficient taken as that, as LAPACK's dlasy2 (the blocks' entries near 1: no scaling against overflow).
|
inlineprivate |
Decomposes the leaf node by implicit-shift QR.
Definition at line 5740 of file decomposition.hpp.
|
inlineprivate |
Decomposes the block index of nodes and everything below it, each block right after its halves.
Definition at line 5885 of file decomposition.hpp.
Solves [a b; c d] x = r in place of r by Gaussian elimination with complete pivoting, a pivot below floor taken as floor, as LAPACK's dlaln2; r is scaled first, and.
x would pass back_substitution_limit.
|
inline |
Solves the leaf node by implicit QL: its eigenvalues ascending into its diagonal, their eigenvectors as the columns of its block of q.
Definition at line 4057 of file decomposition.hpp.
|
inline |
Solves the block index of nodes and everything below it, each block right after its halves.
Definition at line 4112 of file decomposition.hpp.
| 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.
The root lies in (Pᵢ, Pᵢ₊₁), or for the last within 2ρ‖z‖² past its pole. It is sought as an offset τ from the nearer end P_o, so that every Pⱼ − x = offset(o, j) − τ is computed with a small relative error, which keeps the vectors built from them orthogonal. Each step fits the equation with its two nearest poles (Gragg's) and falls back to bisection when the fit leaves the bracket, which shrinks at every step.
| [in] | z | The weights. |
| [in] | rho | ρ. |
| [in] | i | The root. |
| [in] | max_iterations | The bound on the steps. |
| [out] | delta | Pⱼ − x for every j (size z.size()). |
| [in] | offset | The pole offsets. |
| convergence_error | If the root is not reached within max_iterations steps. |
Definition at line 3613 of file decomposition.hpp.
| void numerical::detail::solve_upward | ( | const buffer< T > & | t, |
| std::size_t | n, | ||
| std::size_t | k, | ||
| std::size_t | top, | ||
| typename X::scalar | lambda, | ||
| double | floor, | ||
| X & | x | ||
| ) |
Solves (T − λ I) x = 0 upward from the component k, rows k − 1 … 0 of the row-major quasi-triangular t (its 2 × 2 blocks where a subdiagonal entry is nonzero), x [k … top] given: as LAPACK's dtrevc and ztrevc, a pivot below floor taken as floor, the components scaled down where they would grow past back_substitution_limit. The sums of rows against x are X 's, vectorised.
Polls ready, yielding the core between polls, until it holds or spin_time has passed.
Definition at line 225 of file parallel.hpp.
|
inlinenoexcept |
Whether gemm_split runs a product of m rows of A by four real products of each block's parts (split_parts_block), else by one of its expanded operands (split_expanded_block): the same bits, the parts' products faster where they read B's parts in place (gemm_packed). On the AMX path up to split_parts_rows; elsewhere (the NEON path) up to one row of the NEON tile, its b_in_place_rows: M1 Pro, all threads, by 400 × 500, 1 to 4 rows 2 to 18 % faster, 6 to 32 rows 5 to 52 % slower.
| void numerical::detail::split_expanded_block | ( | std::size_t | m, |
| std::size_t | nb, | ||
| std::size_t | kb, | ||
| std::size_t | pc, | ||
| std::size_t | jc, | ||
| const Z * | a, | ||
| std::size_t | lda, | ||
| const Z * | b, | ||
| std::size_t | ldb, | ||
| typename Z::value_type * | ar, | ||
| typename Z::value_type * | br, | ||
| typename Z::value_type * | cr | ||
| ) |
One block of gemm_split (split_parts_block's) by one real product of its expanded operands, [Re C | Im C] += [Re A | −Im A] · [[Re B, Im B], [Im B, −Re B]]: the same chains, B copied four times its size, the product as wide and as deep as the block's two.
|
inlineprivate |
Splits an einsum left-hand side into per-operand subscripts.
Definition at line 347 of file tensor.hpp.
| void numerical::detail::split_parts_block | ( | std::size_t | m, |
| std::size_t | nb, | ||
| std::size_t | kb, | ||
| std::size_t | pc, | ||
| std::size_t | jc, | ||
| const Z * | a, | ||
| std::size_t | lda, | ||
| const Z * | b, | ||
| std::size_t | ldb, | ||
| typename Z::value_type * | ar, | ||
| typename Z::value_type * | br, | ||
| typename Z::value_type * | cr | ||
| ) |
One block of gemm_split ([Re C | Im C] at cr, m × nb each, += A·B, A m × kb and B kb × nb) by four real products of its parts: Re C += Re A·Re B, then −Im A·Im B; Im C += Re A·Im B, then Im A·Re B, each continuing its elements' chains ((−x)·(−y) is x·y exactly: the bits of −Im A·−Re B). B's parts are a copy of its size, which a product of a few rows of A reads in place. The block is A's columns and B's rows from pc, B's columns from jc, A and B laid out as A and B say.
| split_operands< R > & numerical::detail::split_storage | ( | ) |
The calling thread's split_operands. A product uses them whole on one thread, and the real product it runs on them nests no complex one.
NumPy's complex square, (re·re − im·im) + (re·im + im·re)i, without the infinity handling of a general product (C Annex G), which costs a call per element.
Definition at line 324 of file tensor_ops.hpp.
x · x.
Definition at line 316 of file tensor_ops.hpp.
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.
| convergence_error | If the iteration does not converge within its bound (not observed). |
Definition at line 6409 of file decomposition.hpp.
|
inline |
|
inline |
|
inlineprivate |
Its singular values and vectors into node's, from scratch: u zero, v the identity.
Definition at line 5728 of file decomposition.hpp.
| 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.
Definition at line 2946 of file decomposition.hpp.
| 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.
Definition at line 2964 of file decomposition.hpp.
|
inline |
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).
Definition at line 669 of file decomposition.hpp.
| 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.
Definition at line 655 of file decomposition.hpp.
|
inline |
y[i] −= f · x[i] for double on the widest x86-64 units the matrix product's kernel may use (gemm_path_current), each element one fused multiply-add there; elsewhere subtract_product's template. The decompositions' row updates take it: a wheel built for the baseline x86-64 holds two double per vector and fuses none.
|
inline |
y[i] −= f · x[i] for std::complex<double>: one term of subtract_products on AVX-512 or AVX2 where the matrix product's kernel runs there, else the template.
y[i] -= f · x[i] for count elements of y and x, which do not overlap. Each element is computed exactly as that scalar statement computes it (its product rounded before its difference, or fused with it where the compiler contracts the statement); the arguments follow the expression, since a compiler may contract f · x and x · f of complex numbers differently. Under GCC, over native vectors where T has them: at -O2 GCC leaves the scalar loop scalar, unable to rule out an overlap without a check it will not add. Clang vectorises the scalar loop itself, with that check, and better than these vectors (an LU solve of order 128 on arm64 took 19 % longer through them). For double, gemm.hpp's overload chooses the x86-64 units at run time, a fused multiply-add per element there.
y[i] -= x[i] · f, as the other subtract_product (a real product is the same either way round).
|
inline |
y[i] -= Σₜ f[t] · x[(t · ldx) + i] for double, the units chosen as subtract_product chooses them, the elements held in registers through the terms there.
|
inline |
subtract_products for std::complex<double>: on AVX-512 or AVX2 with FMA where the matrix product's kernel runs there, each part of each term one fused multiply-add (Re f · x, then Im f times the swapped parts), the elements in registers through the terms; elsewhere the template, whole complex products as the compiler rounds them.
| void numerical::detail::subtract_products | ( | T * | y, |
| const T * | f, | ||
| const T * | x, | ||
| std::size_t | ldx, | ||
| std::size_t | terms, | ||
| std::size_t | count | ||
| ) |
y[i] -= Σₜ f[t] · x[(t · ldx) + i] for count elements of y, terms rows of x apart by ldx: subtract_product for each term in ascending order, so each element is computed as that sequence of calls computes it. For double, gemm.hpp's overload holds the elements in registers through the terms.
| void numerical::detail::subtract_products | ( | T * | y, |
| const T * | f, | ||
| const T * | x, | ||
| std::size_t | ldx, | ||
| std::size_t | terms, | ||
| std::size_t | count | ||
| ) |
subtract_products for a real T: eight elements of y held in a local array through the terms, which the compiler keeps in registers, where each term was a call through the whole row. Each element takes its terms in the same order as a sequence of subtract_product gives them (an LU solve of order 64 and four columns on arm64: 18 → 8.5 µs), and the same value under Clang, and under GCC on x86-64, which fuses neither; under GCC on arm64 the vectors of subtract_combined_vectors and subtract_scaled_vectors round apart and the scalar statements fused, and their splits differ.
|
inline |
subtract_rows for double, the units chosen as subtract_products chooses them: each element the value subtract_products gives it.
| void numerical::detail::subtract_rows | ( | T * | y, |
| std::size_t | ldy, | ||
| const T * | f, | ||
| std::size_t | ldf, | ||
| const T * | x, | ||
| std::size_t | ldx, | ||
| std::size_t | terms, | ||
| std::size_t | count | ||
| ) |
subtract_products for product_rows rows of y at once, ldy apart, each with its own factors (ldf apart in f) on the same rows of x: y_r[i] -= Σₜ f_r[t] · x[(t · ldx) + i]. Each element takes the statement of subtract_products, a term at a time in ascending order, so the same value; each element of x read serves every row, and the rows' chains of multiply-adds advance together, where a row alone waits on its last one through each term.
| 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).
Definition at line 99 of file decomposition.hpp.
sum −= a · b in one statement, which the compiler may fuse (as it fuses the matrix solve's).
Definition at line 90 of file decomposition.hpp.
|
private |
einsum of one row-major operand data (shape, its axes labelled in, each label once) into the labels output, where those it sums are one run of its axes and output the others in order (ij->i, ij->, ij->j, ijk->ik): each output element the sum of evenly spaced elements, by pairwise halving as sum() adds them (pairwise_sums), where einsum's general loop added them one after another. None for any other spec (nothing summed, the summed labels apart, the output reordered).
Definition at line 500 of file tensor.hpp.
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).
Definition at line 2851 of file decomposition.hpp.
| std::array< T, symmetric_row_group > numerical::detail::symmetric_rows | ( | const T * | rows, |
| std::size_t | n, | ||
| const T * | v, | ||
| T * | out, | ||
| std::size_t | i | ||
| ) |
symmetric_rows one column at a time, for a type no native vector holds.
| std::array< T, symmetric_row_group > numerical::detail::symmetric_rows_tail | ( | const T * | rows, |
| std::size_t | n, | ||
| const T * | v, | ||
| T * | out, | ||
| std::size_t | i, | ||
| std::size_t | done, | ||
| std::array< T, symmetric_row_group > | sums | ||
| ) |
The columns done…i+r of symmetric_rows, one by one: for each row r in order, out[k] += row_r[k] · v[i+r] for done ≤ k < i + r and row_r[k] · v[k] added to sums[r], its diagonal term last.
|
constexprnoexcept |
task_rows_cap for a tile whose task rows depend on the panel's depth (Kernel::task_rows_for).
|
constexprnoexcept |
|
inlineprivate |
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.
Definition at line 5402 of file decomposition.hpp.
|
inlinenoexcept |
The most threads a parallel_for uses (0: as many as the pool has), set by numerical::set_max_threads.
Definition at line 525 of file parallel.hpp.
|
inline |
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.
Definition at line 3761 of file decomposition.hpp.
|
inline |
How many threads count indices of work_per_index each pay for: one per work_per_thread, at most one per index and per hardware thread.
| [in] | count | The number of indices. |
| [in] | work_per_index | The work of one index, in multiply-adds. |
| [in] | hardware | The hardware threads. |
Definition at line 553 of file parallel.hpp.
|
noexcept |
One tile of a matrix product kernel, rows × cols of C (at most Rows × Cols, leading dimension ldc), through full, which computes a whole tile in place given its first element and leading dimension: a whole tile in C itself, a partial one in a copy whose outside is zero. Only that copy is zeroed: zeroing it for every tile, whole ones included, took a twelfth of an x86-64 product's time.
|
private |
Matrix product C(M×N) = A(M×K) · B(K×N), row-major.
The fast path for a plain two-operand single-contraction einsum (matmul, matrix–vector, dot), by detail::gemm. Every element of C is summed in ascending order of the contraction, the general engine's order, so no tolerance changes.
Definition at line 375 of file tensor.hpp.
|
inline |
a · b of two complex numbers by the schoolbook formula, without C Annex G's recovery of infinities, which a large translation unit leaves out of line (a call per product: complex eig of order 64 took 2.7 times as long in the Python binding as alone). The entries here are finite.
| 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.
Definition at line 3327 of file decomposition.hpp.
|
private |
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).
Definition at line 2235 of file decomposition.hpp.
| 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(δⱼ).
Definition at line 3300 of file decomposition.hpp.
| 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.
| convergence_error | If an iteration did not converge within its bound. |
Definition at line 4239 of file decomposition.hpp.
| 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.
| [in,out] | d | The diagonal; on return, the eigenvalues (unsorted). |
| [in,out] | e | The subdiagonal in e[0…n-2], e[n-1] = 0; destroyed. |
| [in,out] | zt | Row-major n × n whose rows are rotated together (Qᵀ from tridiagonalize); on return, row j is the eigenvector of d[j]. |
| [in] | n | The order. |
| [in] | max_steps | The bound of QL steps in all. |
| convergence_error | If a subdiagonal entry is not negligible after max_steps steps. |
Definition at line 3491 of file decomposition.hpp.
| 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.
Each reflection zeroes one column below the subdiagonal; the trailing block takes the rank-2 updates B − v wᴴ − w vᴴ of a panel of tridiagonal_block reflections at once, through the matrix product kernel, as LAPACK's dsytrd. A complex subdiagonal entry is then made real by a unit phase per row, as LAPACK's zhetrd leaves it for a real iteration; a real one keeps its sign.
| [in,out] | a | The matrix; destroyed. |
| [in] | n | Its order (at least 1). |
| [out] | d | The diagonal of T (size n). |
| [out] | e | The subdiagonal of T in e[0…n-2], and e[n-1] = 0 (size n). |
M's reflectors and phases: tridiagonal_basis forms M, to_matrix_basis applies it. Definition at line 3174 of file decomposition.hpp.
|
inline |
|
inline |
|
inline |
The eigenvalues of the real Schur form t (jw × jw) of a deflation window that do not deflate, moved to its top, as dlaqr3: from the bottom, a block whose spike entries s v₁ⱼ are negligible deflates, any other is moved up; the undeflated blocks are then sorted by decreasing modulus (which helps graded matrices), and every eigenvalue of the window written to w.
|
inline |
|
inlineprivate |
The singular values of a solved node scaled back.
Definition at line 5863 of file decomposition.hpp.
|
inline |
The eigenvalues of a solved node, scaled back, where its parent reads them.
Definition at line 4102 of file decomposition.hpp.
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.
Definition at line 6296 of file decomposition.hpp.
|
private |
canonical_product's vector a (inner) times a matrix given as its transpose bt (cols × inner, its rows ldbt apart, B's columns): each element of the product the sum canonical_product forms, a[k] · B[k][j] added in ascending order of k with add_product's arithmetic (add_one_product), so the same bits; a column of B a run of bt, transposed_columns of them at once.
Definition at line 439 of file tensor.hpp.
| void numerical::detail::vigilant_deflation | ( | one_based< T > | h, |
| lapack_index | k, | ||
| lapack_index | ktop, | ||
| lapack_index | kbot, | ||
| double | smallest | ||
| ) |
|
inline |
axis_walk of shape along axis; outer and inner 0 where the result is empty.
| dimension_error | If axis is not an axis, or the result's elements overflow. |
Definition at line 71 of file tensor_ops.hpp.
|
constexprnoexcept |
Whether the worker on core joins a run of helpers workers whose caller is on caller: by rank among the cores, the caller's aside, so a run of few threads takes the same cores from one run to the next.
Definition at line 291 of file parallel.hpp.
|
inline |
The shape with axis axis removed.
Definition at line 287 of file tensor_ops.hpp.
|
inline |
The core of each worker the pool starts beside its callers, from the cores known (0: none) and its threads (at least 1): one worker per core, pinned there, where the cores are known (each run then leaves out the one on the caller's core, sits_out); else threads − 1 workers free to move.
Definition at line 66 of file parallel.hpp.
|
inline |
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).
Definition at line 3737 of file decomposition.hpp.
|
inlineconstexprprivate |
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.
Definition at line 4755 of file decomposition.hpp.
|
inlineconstexprprivate |
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).
The reflections from each side in one panel of the blocked bidiagonal reduction.
Definition at line 4816 of file decomposition.hpp.
|
inlineconstexprprivate |
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 %.
Definition at line 4750 of file decomposition.hpp.
|
inlineconstexprprivate |
The order up to which bidiagonal_divide decomposes by implicit-shift QR.
Definition at line 5657 of file decomposition.hpp.
|
inlineconstexprprivate |
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 %.
Definition at line 4824 of file decomposition.hpp.
|
inlineconstexpr |
The columns of one panel of the blocked Cholesky factorisation (as lu_block).
Definition at line 200 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 204 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 216 of file decomposition.hpp.
|
inlineconstexpr |
The rows of one chunk of a Cholesky panel's rows below it, solved together along each column.
Definition at line 235 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 208 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 240 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 245 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 2940 of file decomposition.hpp.
|
inlineconstexpr |
The order up to which tridiagonal_divide_and_conquer diagonalises by implicit QL.
Definition at line 3565 of file decomposition.hpp.
|
inlineconstexpr |
The elements one task of an element-wise operation takes, a quarter of a thread's work: a tensor of a million shares out (M1 Pro, 2000 × 2000 doubles added: 3.1 → 1.4 ms), one of a few hundred thousand stays within its time (a sixteenth: 300 × 300 plus a row 13 % slower, each task setting up its odometer).
Definition at line 539 of file parallel.hpp.
|
inlineconstexpr |
|
inlineconstexprprivate |
The operands at most whose contraction order contraction_path searches exhaustively; more take the cheapest pair at each step.
Definition at line 688 of file tensor.hpp.
The elements all_finite_shared gives a thread at a time.
Definition at line 553 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 3729 of file decomposition.hpp.
|
inlineconstexpr |
|
inlineconstexpr |
The multiply-adds the register tile does in the time of one in work_per_thread's dependent chain (about 20 G/s on one arm64 or x86-64 AVX2 core): a task weighs its work by it, so a product of 0.1 ms does not start eight threads.
|
inlineconstexpr |
|
inlineconstexpr |
|
inlineconstexpr |
|
inlineconstexpr |
Below this many multiply-adds the plain loop (gemm_rows, by blocks of rows held in registers) runs instead of the packed kernel, both summing in the same order. Against the packed kernel on 400 random shapes from 512 multiply-adds, the plain loop took 0.25 to 0.41 times as long below 2^13 on arm64 (M1 Pro, geometric means by octave, the slowest shape 0.91); above, it lost up to 2.5 times on shapes the AMX suits. On x86-64 the plain loop has the baseline's vectors and the packed one the AVX2 or AVX-512 units: there it took 0.43 to 0.55 times as long below 2^11, 0.86 above with shapes up to 2 times slower (six AVX2 cores, GCC 13).
|
inlineconstexpr |
|
inlineconstexpr |
The columns of a panel of hessenberg_blocked.
|
inlineconstexpr |
The active order from which hessenberg reduces by panels.
|
inlineconstexpr |
The rows an LU inverse writes its identity into, and permutes the columns of, on one thread at a time.
Definition at line 174 of file decomposition.hpp.
The |ζ| past which one-sided Jacobi takes the rotation's tangent as 1 / (2|ζ|): ζ² stays finite under it.
Definition at line 4635 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 2687 of file decomposition.hpp.
|
inlineconstexprprivate |
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 / ε).
Definition at line 4632 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 3754 of file decomposition.hpp.
|
inlineconstexpr |
The poles a task of merge_rank_one's Löwner products takes, and the side of its tiles of u.
Definition at line 3708 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 108 of file decomposition.hpp.
|
inlineconstexpr |
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.
Definition at line 136 of file decomposition.hpp.
|
inlineconstexpr |
The columns of one panel of the LU of an order up to lu_small_order.
Definition at line 124 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 120 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 144 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 178 of file decomposition.hpp.
|
inlineconstexpr |
|
inlineconstexpr |
|
inlineconstexpr |
|
inlineconstexpr |
The rows of C one block of gemm_rows takes together: each element a chain of multiply-adds that waits on its last, four rows give the units four times the chains (M1 Pro, order 8 to 24: 1.6 to 1.9 times as fast as one row; two rows on six x86-64 AVX2 cores, the fastest there, were within 7 %).
|
inlineconstexpr |
The least depth at which rows of C wider than one block (narrow_columns) run by blocks in gemm_rows; one block or narrower takes any depth past 1. Against the plain loop of unknown width, on 500 random shapes below the packed kernel, the blocks took 0.11 to 0.63 times as long from a depth of 2 on arm64 (M1 Pro, geometric means by width); on x86-64, where the loop has the baseline's vectors, 0.35 to 0.90 from a depth of 8 but 1.08 to 1.31 below it past one block (six AVX2 cores, GCC 13).
|
inlineconstexpr |
The most right-hand sides an LU solve takes row by row rather than by blocks through the product kernel.
Definition at line 148 of file decomposition.hpp.
|
inlineconstexpr |
|
inlineconstexpr |
No core: a worker free to run anywhere, or a caller on a processor the pool does not know.
Definition at line 59 of file parallel.hpp.
|
inlineconstexpr |
|
inlineconstexprprivate |
The columns pairwise_sum_rows sums at once: its partial sums of eight rows, 32 KiB of doubles, stay in the L1 cache.
Definition at line 314 of file tensor.hpp.
|
inlineconstexpr |
The fewest columns of a product plain_loop_pays sends to the plain loop on the x86-64 tiles: there an order-24 SVD, whose two products of 24 × 24 × 24 it took, ran 9 % slower (GCC 13, the same bits), though such a product alone ran faster in the plain loop. On the NEON tile the small square products gain most (GCC 14, one core: a product of order 24 to 64 −32 to −50 %, a QR −5 to −21 %).
The largest plain sum of squares plain_or_at_scale keeps: no term of it overflowed.
Definition at line 373 of file decomposition.hpp.
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.
Definition at line 370 of file decomposition.hpp.
Whether Kernel has the C tile of the next row of tiles prefetched while a tile computes (Kernel::prefetch_rows): the tiles within a row walk one page's lines, which the processor's own prefetch follows, while the next row lands on other pages.
|
inlineconstexpr |
The rows of y subtract_rows takes at once.
|
inlineconstexpr |
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.
Definition at line 3123 of file decomposition.hpp.
|
inlineconstexprprivate |
The reflectors of one panel of the blocked QR.
Definition at line 2143 of file decomposition.hpp.
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.
Definition at line 6271 of file decomposition.hpp.
|
inlineconstexprprivate |
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.
Definition at line 2230 of file decomposition.hpp.
|
inlineconstexpr |
The magnitude under which rescue_factor scales a T up: about the square root of the smallest normal number of its range.
Definition at line 301 of file decomposition.hpp.
|
inlineconstexpr |
The rows a thread takes of a matrix times a vector (matrix_times_vector).
|
inlineconstexpr |
|
inlineconstexpr |
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.
Definition at line 3569 of file decomposition.hpp.
|
inlineconstexpr |
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).
Definition at line 159 of file decomposition.hpp.
|
inlineconstexpr |
The columns of one leaf of the recursive LU panel.
Definition at line 171 of file decomposition.hpp.
|
inlineconstexpr |
How long a thread polls before it sleeps: kernels run in quick succession (a matrix product packs, then computes, panel after panel), and a thread woken from sleep starts tens of microseconds late, its wake-up paid by the caller of the next run, about 10 µs a sleeping worker on six x86-64 cores. A factorisation's steps last about 100 µs at order 256, so at 100 µs the workers one step left out were asleep for the next (an LU of order 256: six x86-64 cores 560 → 500 µs at 500; eight arm64 cores 358 → 324).
Definition at line 221 of file parallel.hpp.
|
inlineconstexpr |
|
inlineconstexpr |
The depth of one real product of gemm_split (split_columns).
|
inlineconstexpr |
|
inlineconstexpr |
The most rows of A for which gemm_split runs four real products of a block's parts on the AMX path (split_by_parts). M1 Pro, all threads, A of 16 to 64 rows by B of 256 to 1024: complex float 1.1 to 2.4 times faster, complex double 0.97 to 1.6; 96 to 128 rows by 512, complex double 1.2 to 1.4 times slower, where one wide product suited the AMX.
|
inlineconstexpr |
The rows of an SVD's input one thread scales and copies at a time.
Definition at line 181 of file decomposition.hpp.
|
inlineconstexprprivate |
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.
Definition at line 4745 of file decomposition.hpp.
|
inlineconstexpr |
The rows of a lower triangle symmetric_rows takes at once.
|
inlineconstexpr |
|
inlineconstexpr |
The side of the tiles a transposing copy goes by (an einsum permutation of two axes, the eigenvectors' columns): a tile read and a tile written fit the L1 cache, where element by element each write, a row apart, missed it (order 1024 on arm64: 4.7 → 1.5 ms; the eigenvectors of order 1024, a tenth of an eigh on six x86-64 cores).
|
inlineconstexprprivate |
The columns of B one task of vector_times_transposed sums together: their sums' chains interleaved, each waiting on its own last addition.
Definition at line 430 of file tensor.hpp.
|
inlineconstexpr |
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.
Definition at line 3119 of file decomposition.hpp.
|
inlineconstexpr |
The reflections of one panel of the blocked tridiagonal reduction.
Definition at line 3114 of file decomposition.hpp.
|
inlineconstexprprivate |
The columns a thread takes of a vector times a matrix: each of its rows added to them in order.
Definition at line 391 of file tensor.hpp.
|
inlineconstexpr |
The multiply-adds a vectorised inner loop does in the time of one in work_per_thread's dependent chain: a product weighs its work by it (at order 64, a matrix product on eight threads took 95 µs, on one 35 µs).
Definition at line 544 of file parallel.hpp.
|
inlineconstexpr |
|
inlineconstexpr |
The work one thread must have to pay for starting it: about 60 µs of multiply-adds, six times the cost of a start (measured on arm64: 16 indices of 4096 ran 80 µs sequentially, 92 µs on eight threads; of 16384, 324 µs and 128 µs).
Definition at line 534 of file parallel.hpp.
|
inlineconstexpr |
The elements one index of parallel_zero writes: eight 4 KiB pages of double.
Definition at line 644 of file parallel.hpp.
|
inlineconstexpr |
The work of zeroing one element, in multiply-adds: a second thread from 1 MiB of double. Below, the allocator mostly hands back memory already touched, and zeroing 512 KiB took 4.5 µs on one arm64 thread, 9 µs split across eight.
Definition at line 671 of file parallel.hpp.
|
inlineconstexpr |
The rows of a Cholesky factor whose upper part one thread zeroes at a time.
Definition at line 184 of file decomposition.hpp.