|
NumericAL
A header-only C++20 lexer built on REAL
|
An LU factorisation with partial pivoting of a square matrix. More...
#include <decomposition.hpp>
Public Member Functions | |
| lu_decomposition (const matrix< T > &a) | |
Factorises a (in place, into combined L\U with a row permutation). | |
| lu_decomposition (matrix< T > &&a) | |
Factorises a in its own storage, an expiring matrix's elements taken over: no copy of them. | |
| vector< T > | solve (const vector< T > &b) const |
Solves A x = b using the stored factors. | |
| T | determinant () const |
| The determinant (the pivot sign times the product of U's diagonal; 0 if singular). | |
| std::pair< T, double > | log_determinant () const |
The determinant as its unit-modulus sign and the natural logarithm of its magnitude, as numpy.linalg.slogdet gives them: a sum of logarithms does not underflow or overflow where the product does (det(1e-2·I) of order 200 is 0.0; its logarithm is −921.0). | |
| matrix< T > | inverse () const |
The inverse, U⁻¹ L⁻¹ P, by blocks of rows through the matrix product kernel: L⁻¹ from the identity, whose row i stays zero past column i in the forward pass, so that pass reads only the triangle it fills (a third of the products of a full one, as LAPACK's dgetri spares them); then U⁻¹ on it, and P as a permutation of the columns. The terms left out are products with zeros. | |
| matrix< T > | inverse (buffer< T > storage) const |
inverse in storage, its elements discarded: no fresh n × n mapped and faulted in where it holds that much already. pivoted_inverse, then P applied in place, by bands of rows on the threads. | |
| matrix< T > | pivoted_inverse (buffer< T > storage) const |
(P A)⁻¹ = U⁻¹ L⁻¹, the inverse of the rows in their pivoted order, in storage (its elements discarded): column pivots()[i] of A⁻¹ is its column i, which inverse moves there (a caller that copies the result anyway can move them as it copies). The identity written by bands of rows on the threads. | |
| matrix< T > | solve (const matrix< T > &b) const |
Solves A X = B for every column of b at once. | |
| matrix< T > | solve (matrix< T > &&b) const |
solve on an expiring b, in its own storage: its rows permuted in place (cycle by cycle of the pivots, one row held aside), then solved there. No copy of B. | |
| matrix< T > | factors () const & |
L and U in one matrix, as LAPACK's dgetrf leaves them: U on and above the diagonal, L below it (its unit diagonal not stored), of the rows in their pivoted order. | |
| matrix< T > | factors () && |
| factors, moved out of an expiring decomposition: no copy. | |
| const std::vector< std::size_t > & | pivots () const noexcept |
The row permutation: row i of P A = L U is row pivots()[i] of A. | |
| bool | singular () const noexcept |
Whether a pivot was zero: U is singular, and solve and inverse throw. | |
Private Member Functions | |
| void | solve_rows (buffer< T > &y, std::size_t cols, bool lower=false) const |
| Performs the pivoted elimination, filling lu_ and pivots_. | |
| std::size_t | column_slices (std::size_t cols) const |
The column slices solve_rows splits cols right-hand sides into, one a thread, each of detail::slice_columns or more (on the x86-64 units, to detail::slice_solve_order): the blocks of rows then join the threads once, where each block joined them for products too small to share (six AVX2 cores, an inverse of order 256: 1.87 → 1.60 ms; arm64, eight cores, 698 → 503 µs). | |
| std::size_t | lower_slice_start (std::size_t s, std::size_t slices, std::size_t cols) const |
The first column of slice s of slices (slices itself: cols) where the forward pass is lower: column c weighs n² + (n − c)², the backward pass's triangle and the forward pass's below row c, and the slices weigh alike (by columns alike the first slice did twice the last's work), each from a multiple of 8 columns. A slice is never empty: it spans at least half of cols / slices columns, past 8. | |
| void | solve_columns (buffer< T > &y, std::size_t cols, std::size_t c0, std::size_t c1, bool lower) const |
solve_rows on the columns [c0, c1) of y (n × cols, row-major): the blocked passes, their products through the matrix kernel, each element's terms in the same order whatever the slice. | |
| void | solve_narrow (buffer< T > &y, std::size_t cols) const |
| solve_rows for a few right-hand sides: row by row, each one's terms in the order the blocked passes give them (forward, ascending; backward, the blocks from the last down to its own, each ascending), without a product that packs the factor's blocks for a few columns (order 256 and four columns on one x86-64 core: 112 → 72 µs, the factorisation 750). The solution is the blocked passes' to the bit wherever their products run the matrix kernel; one the plain loop takes (under detail::gemm_packed_threshold multiply-adds, or shallower than detail::gemm_min_depth, as a last block of a few rows) rounds each term apart where the x86-64 rows fuse it, so there the two differ in the last bits. | |
| std::size_t | panel_of (std::size_t i) const noexcept |
The first row of the panel that holds row i. | |
| void | subtract_block (buffer< T > &y, std::size_t cols, std::size_t r0, std::size_t r1, std::size_t c0, std::size_t c1, std::size_t first, std::size_t last) const |
Y[r0, r1) −= F[r0, r1) × [c0, c1) · Y[c0, c1), F the factor, through the matrix product kernel (nothing when the rows are empty), on the columns [first, last) of Y, first < last. | |
| void | negated_block_into (buffer< T > &out, std::size_t r0, std::size_t r1, std::size_t c0, std::size_t c1) const |
negated_block into out, resized to it. | |
| buffer< T > | negated_block (std::size_t r0, std::size_t r1, std::size_t c0, std::size_t c1) const |
The rows [r0, r1) and columns [c0, c1) of the factor, negated, row-major and contiguous. | |
| void | factorize () |
| void | factorize_trailing (buffer< T > &panel, std::size_t first_end) |
The panels after the first (whose columns end at first_end), each with its trailing update. Out of line: in factorize, the update's storage and calls cost an LU of one panel 5 % (orders 8 to 48 under clang on arm64), where none of them runs. | |
| void | factorize_panel (buffer< T > &panel, std::size_t k0, std::size_t end) |
Unblocked LU with partial pivoting of the columns [k0, end), every row below k0. | |
| void | transpose_panel (buffer< T > &panel, std::size_t rows, std::size_t k0, std::size_t end, bool in) |
Copies the columns [k0, end) of the factor, every row from k0, into panel (column-major, rows rows) when in, back otherwise: by tiles of rows on the x86-64 units from panel_tile_rows rows, row by row otherwise, in the plainest loops, which Apple clang compiles best (tiles of one row there ran the whole LU of order 17 1.1 → 1.7 µs, of order 512 1.8 → 2.1 ms, on one thread). | |
| void | factorize_columns (buffer< T > &panel, std::size_t rows, std::size_t first, std::size_t last, std::size_t k0, std::size_t end) |
Factorises the columns [first, last) of the panel (column-major, rows rows), every row from first down, each pivot's row swap applied to the whole panel and outside it. | |
| void | swap_rows_outside (std::size_t a, std::size_t b, std::size_t k0, std::size_t end) |
Swaps rows a and b of the factor outside the panel's columns [k0, end): right of the panel up to swap_limit_, the rest left in deferred_swaps_ (columns another thread is updating); left of it, in left_swaps_, for swap_left_columns. | |
| void | swap_left_columns () |
| Applies left_swaps_ to the columns left of each swap's panel: a panel's columns take, in order, the swaps of the panels after it, which are the swaps from its last column on. The panels' columns are independent, so they go to the threads; a column sees its swaps in the order the pivots were chosen. | |
| void | solve_panel_rows (std::size_t k0, std::size_t end) |
U₁₂ = L₁₁⁻¹ A₁₂: the panel's rows right of it, by forward substitution with its unit-lower L. Each block of columns is independent, so blocks run in parallel, each in the sequential order: each element still takes its terms in ascending order of the pivot. | |
| void | update_columns (const detail::packed_left< T > &left, std::size_t k0, std::size_t c0, std::size_t c1) |
A₂₂ −= L₂₁ U₁₂ on the trailing columns [c0, c1), as A₂₂ += (−L₂₁) U₁₂ through the matrix product kernel, left being −L₂₁ (the rows below the panel from k0) as detail::pack_left packed it. | |
Private Attributes | |
| std::size_t | n_ |
| The matrix order. | |
| std::size_t | block_ {detail::lu_block_for(n_)} |
| The panel width. | |
| buffer< T > | lu_ |
| Combined L (below diag, unit) and U (on/above), row-major. | |
| std::vector< std::size_t > | pivots_ |
| Row permutation: row now at i came from pivots_[i]. | |
| int | sign_ {1} |
| Permutation sign, for the determinant. | |
| bool | singular_ {false} |
| Set if a zero pivot was met. | |
| std::size_t | swap_limit_ {n_} |
| The columns a panel's row swaps reach now. | |
| std::vector< std::pair< std::size_t, std::size_t > > | deferred_swaps_ |
| Swaps left for the columns past it. | |
| std::vector< std::pair< std::size_t, std::size_t > > | left_swaps_ |
| Swaps left for the columns before their panel. | |
An LU factorisation with partial pivoting of a square matrix.
| T | The scalar element type (a field: double, std::complex<double>, …). |
Factorises once on construction; solve, determinant and inverse then reuse the factors. Built for the common case of one matrix and several right-hand sides.
Definition at line 688 of file decomposition.hpp.
|
inlineexplicit |
Factorises a (in place, into combined L\U with a row permutation).
| [in] | a | The square matrix to factorise. |
| dimension_error | If a is not square. |
Definition at line 697 of file decomposition.hpp.
|
inlineexplicit |
Factorises a in its own storage, an expiring matrix's elements taken over: no copy of them.
| [in] | a | The square matrix to factorise. |
| dimension_error | If a is not square. |
Definition at line 706 of file decomposition.hpp.
|
inlineprivate |
The column slices solve_rows splits cols right-hand sides into, one a thread, each of detail::slice_columns or more (on the x86-64 units, to detail::slice_solve_order): the blocks of rows then join the threads once, where each block joined them for products too small to share (six AVX2 cores, an inverse of order 256: 1.87 → 1.60 ms; arm64, eight cores, 698 → 503 µs).
Definition at line 994 of file decomposition.hpp.
|
inline |
The determinant (the pivot sign times the product of U's diagonal; 0 if singular).
Definition at line 755 of file decomposition.hpp.
|
inlineprivate |
Definition at line 1207 of file decomposition.hpp.
|
inlineprivate |
Factorises the columns [first, last) of the panel (column-major, rows rows), every row from first down, each pivot's row swap applied to the whole panel and outside it.
Recursively, as LAPACK's dgetrf2: the left half is factorised, the top of the right half solved against its unit-lower L, the bottom updated by a product through the matrix product kernel (in parallel, where the work is), then the right half factorised; a leaf of detail::solve_leaf columns, unblocked. Each element still receives its updates in ascending order of the pivot.
Definition at line 1349 of file decomposition.hpp.
|
inlineprivate |
Unblocked LU with partial pivoting of the columns [k0, end), every row below k0.
The panel is factorised in a column-major copy, where a pivot search and every update run along a column at unit stride (row-major, each reads one element per row: four times slower at order 2048); a pivot's row swap reaches the columns outside the panel in place. The arithmetic is the in-place algorithm's, step for step. panel is the copy's storage, of (n − k0) × (end − k0) elements or more.
Definition at line 1283 of file decomposition.hpp.
|
inlineprivate |
The panels after the first (whose columns end at first_end), each with its trailing update. Out of line: in factorize, the update's storage and calls cost an LU of one panel 5 % (orders 8 to 48 under clang on arm64), where none of them runs.
Definition at line 1226 of file decomposition.hpp.
factors, moved out of an expiring decomposition: no copy.
Definition at line 927 of file decomposition.hpp.
L and U in one matrix, as LAPACK's dgetrf leaves them: U on and above the diagonal, L below it (its unit diagonal not stored), of the rows in their pivoted order.
Definition at line 921 of file decomposition.hpp.
The inverse, U⁻¹ L⁻¹ P, by blocks of rows through the matrix product kernel: L⁻¹ from the identity, whose row i stays zero past column i in the forward pass, so that pass reads only the triangle it fills (a third of the products of a full one, as LAPACK's dgetri spares them); then U⁻¹ on it, and P as a permutation of the columns. The terms left out are products with zeros.
| singular_matrix_error | If the matrix is singular. |
Definition at line 794 of file decomposition.hpp.
inverse in storage, its elements discarded: no fresh n × n mapped and faulted in where it holds that much already. pivoted_inverse, then P applied in place, by bands of rows on the threads.
| [in] | storage | The storage of the result, taken over. |
storage. | singular_matrix_error | If the matrix is singular. |
Definition at line 806 of file decomposition.hpp.
|
inline |
The determinant as its unit-modulus sign and the natural logarithm of its magnitude, as numpy.linalg.slogdet gives them: a sum of logarithms does not underflow or overflow where the product does (det(1e-2·I) of order 200 is 0.0; its logarithm is −921.0).
+1 or −1 for a real matrix, a unit complex number for a complex one) and the logarithm; a singular matrix gives a sign of 0 and a logarithm of −∞. Definition at line 771 of file decomposition.hpp.
|
inlineprivate |
The first column of slice s of slices (slices itself: cols) where the forward pass is lower: column c weighs n² + (n − c)², the backward pass's triangle and the forward pass's below row c, and the slices weigh alike (by columns alike the first slice did twice the last's work), each from a multiple of 8 columns. A slice is never empty: it spans at least half of cols / slices columns, past 8.
Definition at line 1013 of file decomposition.hpp.
|
inlineprivate |
The rows [r0, r1) and columns [c0, c1) of the factor, negated, row-major and contiguous.
Definition at line 1178 of file decomposition.hpp.
|
inlineprivate |
negated_block into out, resized to it.
Definition at line 1163 of file decomposition.hpp.
|
inlineprivatenoexcept |
The first row of the panel that holds row i.
Definition at line 1139 of file decomposition.hpp.
(P A)⁻¹ = U⁻¹ L⁻¹, the inverse of the rows in their pivoted order, in storage (its elements discarded): column pivots()[i] of A⁻¹ is its column i, which inverse moves there (a caller that copies the result anyway can move them as it copies). The identity written by bands of rows on the threads.
| [in] | storage | The storage of the result, taken over. |
(P A)⁻¹, in storage. | singular_matrix_error | If the matrix is singular. |
Definition at line 833 of file decomposition.hpp.
|
inlinenoexcept |
The row permutation: row i of P A = L U is row pivots()[i] of A.
Definition at line 933 of file decomposition.hpp.
|
inlinenoexcept |
Whether a pivot was zero: U is singular, and solve and inverse throw.
Definition at line 939 of file decomposition.hpp.
Solves A X = B for every column of b at once.
By blocks of rows, the products of a block's solution into the rows still to solve go through the matrix product kernel. The forward substitution is the vector solve's, term for term; the backward one gives each row the farther blocks' terms first, so a column agrees with the vector solve to rounding.
| [in] | b | The right-hand sides (n × k). |
X (n × k). | dimension_error | If b's row count differs from the matrix. |
| singular_matrix_error | If the matrix is singular. |
Definition at line 862 of file decomposition.hpp.
Solves A x = b using the stored factors.
| [in] | b | The right-hand side. |
x. | dimension_error | If b's size differs from the matrix. |
| singular_matrix_error | If the matrix is singular. |
Definition at line 728 of file decomposition.hpp.
solve on an expiring b, in its own storage: its rows permuted in place (cycle by cycle of the pivots, one row held aside), then solved there. No copy of B.
| [in] | b | The right-hand sides (n × k), taken over. |
X (n × k), in b's storage. | dimension_error | If b's row count differs from the matrix. |
| singular_matrix_error | If the matrix is singular. |
Definition at line 889 of file decomposition.hpp.
|
inlineprivate |
solve_rows on the columns [c0, c1) of y (n × cols, row-major): the blocked passes, their products through the matrix kernel, each element's terms in the same order whatever the slice.
Definition at line 1038 of file decomposition.hpp.
|
inlineprivate |
solve_rows for a few right-hand sides: row by row, each one's terms in the order the blocked passes give them (forward, ascending; backward, the blocks from the last down to its own, each ascending), without a product that packs the factor's blocks for a few columns (order 256 and four columns on one x86-64 core: 112 → 72 µs, the factorisation 750). The solution is the blocked passes' to the bit wherever their products run the matrix kernel; one the plain loop takes (under detail::gemm_packed_threshold multiply-adds, or shallower than detail::gemm_min_depth, as a last block of a few rows) rounds each term apart where the x86-64 rows fuse it, so there the two differ in the last bits.
Definition at line 1088 of file decomposition.hpp.
|
inlineprivate |
U₁₂ = L₁₁⁻¹ A₁₂: the panel's rows right of it, by forward substitution with its unit-lower L. Each block of columns is independent, so blocks run in parallel, each in the sequential order: each element still takes its terms in ascending order of the pivot.
Definition at line 1464 of file decomposition.hpp.
|
inlineprivate |
Performs the pivoted elimination, filling lu_ and pivots_.
L U X = Y in place, y being n × cols, row-major, already permuted. The rows go by blocks of block_ rows: within a block the rows are solved one after another, in parallel chunks of columns (each chunk's columns depend on no other's); from a block to the rows beyond it, through the matrix product kernel. When lower, row i of y is zero past column i (the identity's rows): the forward pass then works on those columns alone. Many columns go to the threads by slices (column_slices), each through every block (solve_columns); a column's terms come in the same order either way.
Definition at line 956 of file decomposition.hpp.
|
inlineprivate |
Y[r0, r1) −= F[r0, r1) × [c0, c1) · Y[c0, c1), F the factor, through the matrix product kernel (nothing when the rows are empty), on the columns [first, last) of Y, first < last.
Definition at line 1146 of file decomposition.hpp.
|
inlineprivate |
Applies left_swaps_ to the columns left of each swap's panel: a panel's columns take, in order, the swaps of the panels after it, which are the swaps from its last column on. The panels' columns are independent, so they go to the threads; a column sees its swaps in the order the pivots were chosen.
Definition at line 1442 of file decomposition.hpp.
|
inlineprivate |
Swaps rows a and b of the factor outside the panel's columns [k0, end): right of the panel up to swap_limit_, the rest left in deferred_swaps_ (columns another thread is updating); left of it, in left_swaps_, for swap_left_columns.
Definition at line 1421 of file decomposition.hpp.
|
inlineprivate |
Copies the columns [k0, end) of the factor, every row from k0, into panel (column-major, rows rows) when in, back otherwise: by tiles of rows on the x86-64 units from panel_tile_rows rows, row by row otherwise, in the plainest loops, which Apple clang compiles best (tiles of one row there ran the whole LU of order 17 1.1 → 1.7 µs, of order 512 1.8 → 2.1 ms, on one thread).
Definition at line 1298 of file decomposition.hpp.
|
inlineprivate |
A₂₂ −= L₂₁ U₁₂ on the trailing columns [c0, c1), as A₂₂ += (−L₂₁) U₁₂ through the matrix product kernel, left being −L₂₁ (the rows below the panel from k0) as detail::pack_left packed it.
Definition at line 1480 of file decomposition.hpp.
|
private |
The panel width.
Definition at line 1490 of file decomposition.hpp.
|
private |
Swaps left for the columns past it.
Definition at line 1496 of file decomposition.hpp.
|
private |
Swaps left for the columns before their panel.
Definition at line 1497 of file decomposition.hpp.
Combined L (below diag, unit) and U (on/above), row-major.
Definition at line 1491 of file decomposition.hpp.
|
private |
The matrix order.
Definition at line 1489 of file decomposition.hpp.
|
private |
Row permutation: row now at i came from pivots_[i].
Definition at line 1492 of file decomposition.hpp.
|
private |
Permutation sign, for the determinant.
Definition at line 1493 of file decomposition.hpp.
Set if a zero pivot was met.
Definition at line 1494 of file decomposition.hpp.
|
private |
The columns a panel's row swaps reach now.
Definition at line 1495 of file decomposition.hpp.
1.9.8