|
NumericAL
A header-only C++20 lexer built on REAL
|
A Cholesky factorisation A = L Lᴴ of a Hermitian positive-definite matrix.
More...
#include <decomposition.hpp>
Public Member Functions | |
| cholesky_decomposition (const matrix< T > &a) | |
Factorises a. | |
| cholesky_decomposition (matrix< T > &&a) | |
Factorises a in its own storage, an expiring matrix's elements taken over: no copy of them. | |
| matrix< T > | lower () const & |
The lower-triangular factor L (with L Lᴴ == A). | |
| matrix< T > | lower () && |
The lower-triangular factor L, moved out of an expiring decomposition. | |
| vector< T > | solve (const vector< T > &b) const |
Solves A x = b via L y = b then Lᴴ x = y. | |
| T | determinant () const |
The determinant — the product of the squared diagonal of L (real, positive). | |
Private Member Functions | |
| void | factorize () |
| void | factor_rows_together (std::size_t i, std::size_t k0) |
The columns [k0, i) of the detail::cholesky_diagonal_rows rows from i, together: each element depends on its own row's earlier columns and on the rows above i alone, so the rows' sums are chains of their own, where row by row each term waited on the last. Each element is factor_row's, its terms in the same order. The rows' loops unrolled: rolled, GCC kept the sums in memory (order 33 on one arm64 core: 8.3 µs, against 2.8 unrolled and 4.6 row by row). | |
| void | factor_row (std::size_t i, std::size_t k0, std::size_t first, std::size_t last) |
Row i of L in the columns [first, last), the terms of the columns before k0 already subtracted: the unblocked step, L[i][j] = (A[i][j] − Σₖ L[i][k]·conj(L[j][k])) / L[j][j] over k ∈ [k0, j). | |
| void | factor_rows_below (std::size_t k0, std::size_t end, buffer< T > &minus_adjoint) |
The rows below the panel [k0, end), L₂₁ = A₂₁ L₁₁⁻ᴴ: chunks of detail::cholesky_rows rows in parallel, each in a column-major copy where column j takes L[i][j] −= L[i][k]·conj(L[j][k]) for k ascending, then the division by L[j][j], for every row of the chunk at once along the column: by sub-blocks of 16 columns, the terms of the sub-blocks before one through the matrix product's plain loop, its own column by column (detail::subtract_product), each element's terms in ascending order, each fused as the column step fuses it, so the unblocked step's arithmetic, in its order (one core of six AVX2, order 1000: 4.95 → 3.95 ms; all column by column, the solve took 37 % of the factorisation for a tenth of its work); row by row, each sum waited on its last addition. Each chunk also writes its columns of minus_adjoint (−L₂₁ᴴ, row-major), the trailing update's right operand, which a copy on one thread after it took 15 to 20 % of the factorisation on six x86-64 cores (order 1024 to 2048). | |
| void | update_trailing (std::size_t k0, std::size_t end, const buffer< T > &minus_adjoint, buffer< T > &packing) |
A22 += L21 (−L21ᴴ), by tiles of detail::cholesky_update_columns columns and detail::cholesky_update_rows rows on and below the diagonal blocks (the upper triangle, never read, takes what falls in a block past the diagonal), in parallel, each element's terms in the order of one product, L₂₁ read in the factor and minus_adjoint as factor_rows_below wrote it: the sign on the right, (−a) b being a (−b) exactly, so no negated copy of L₂₁ (one AVX2 core, order 100: 48 → 45 µs, 300 −4 %). L₂₁ is packed once, in packing, for the tile the products run on (detail::pack_left), and every tile reads its rows there (a tile packed its own, again for each block of columns: 8 % of the factorisation on one AVX2 core at order 1000): so a tile's rows start at a multiple of the tile's slivers, a few rows above the diagonal block. | |
Private Attributes | |
| std::size_t | n_ |
| The matrix order. | |
| std::size_t | block_ {detail::cholesky_block_for(n_)} |
| The panel width. | |
| buffer< T > | lower_ |
| The lower-triangular factor L, row-major (upper triangle zero). | |
A Cholesky factorisation A = L Lᴴ of a Hermitian positive-definite matrix.
| T | The scalar element type (a field: double, std::complex<double>, …). |
L is lower-triangular with a real, positive diagonal. Only the lower triangle of a is read, so a is assumed Hermitian (symmetric, for a real scalar); the factorisation is about twice as cheap as an LU for this special case and cannot be singular. A matrix that is not positive-definite throws not_positive_definite_error.
Definition at line 1584 of file decomposition.hpp.
|
inlineexplicit |
Factorises a.
| [in] | a | The Hermitian positive-definite matrix, its lower triangle alone read. |
| dimension_error | If a is not square. |
| non_finite_error | If an entry of its lower triangle is infinite or NaN. |
| not_positive_definite_error | If a is not positive-definite. |
Definition at line 1594 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 Hermitian positive-definite matrix, its lower triangle alone read. |
| dimension_error | If a is not square. |
| non_finite_error | If an entry of its lower triangle is infinite or NaN. |
| not_positive_definite_error | If a is not positive-definite. |
Definition at line 1605 of file decomposition.hpp.
|
inline |
The determinant — the product of the squared diagonal of L (real, positive).
Definition at line 1659 of file decomposition.hpp.
|
inlineprivate |
Row i of L in the columns [first, last), the terms of the columns before k0 already subtracted: the unblocked step, L[i][j] = (A[i][j] − Σₖ L[i][k]·conj(L[j][k])) / L[j][j] over k ∈ [k0, j).
| not_positive_definite_error | At a diagonal that is not positive. |
Definition at line 1742 of file decomposition.hpp.
|
inlineprivate |
The rows below the panel [k0, end), L₂₁ = A₂₁ L₁₁⁻ᴴ: chunks of detail::cholesky_rows rows in parallel, each in a column-major copy where column j takes L[i][j] −= L[i][k]·conj(L[j][k]) for k ascending, then the division by L[j][j], for every row of the chunk at once along the column: by sub-blocks of 16 columns, the terms of the sub-blocks before one through the matrix product's plain loop, its own column by column (detail::subtract_product), each element's terms in ascending order, each fused as the column step fuses it, so the unblocked step's arithmetic, in its order (one core of six AVX2, order 1000: 4.95 → 3.95 ms; all column by column, the solve took 37 % of the factorisation for a tenth of its work); row by row, each sum waited on its last addition. Each chunk also writes its columns of minus_adjoint (−L₂₁ᴴ, row-major), the trailing update's right operand, which a copy on one thread after it took 15 to 20 % of the factorisation on six x86-64 cores (order 1024 to 2048).
Definition at line 1777 of file decomposition.hpp.
|
inlineprivate |
The columns [k0, i) of the detail::cholesky_diagonal_rows rows from i, together: each element depends on its own row's earlier columns and on the rows above i alone, so the rows' sums are chains of their own, where row by row each term waited on the last. Each element is factor_row's, its terms in the same order. The rows' loops unrolled: rolled, GCC kept the sums in memory (order 33 on one arm64 core: 8.3 µs, against 2.8 unrolled and 4.6 row by row).
Definition at line 1713 of file decomposition.hpp.
|
inlineprivate |
Definition at line 1674 of file decomposition.hpp.
The lower-triangular factor L, moved out of an expiring decomposition.
Definition at line 1623 of file decomposition.hpp.
The lower-triangular factor L (with L Lᴴ == A).
Definition at line 1617 of file decomposition.hpp.
Solves A x = b via L y = b then Lᴴ x = y.
| [in] | b | The right-hand side. |
x. | dimension_error | If b's size differs from the matrix. |
Definition at line 1634 of file decomposition.hpp.
|
inlineprivate |
A22 += L21 (−L21ᴴ), by tiles of detail::cholesky_update_columns columns and detail::cholesky_update_rows rows on and below the diagonal blocks (the upper triangle, never read, takes what falls in a block past the diagonal), in parallel, each element's terms in the order of one product, L₂₁ read in the factor and minus_adjoint as factor_rows_below wrote it: the sign on the right, (−a) b being a (−b) exactly, so no negated copy of L₂₁ (one AVX2 core, order 100: 48 → 45 µs, 300 −4 %). L₂₁ is packed once, in packing, for the tile the products run on (detail::pack_left), and every tile reads its rows there (a tile packed its own, again for each block of columns: 8 % of the factorisation on one AVX2 core at order 1000): so a tile's rows start at a multiple of the tile's slivers, a few rows above the diagonal block.
Definition at line 1857 of file decomposition.hpp.
|
private |
The panel width.
Definition at line 1901 of file decomposition.hpp.
The lower-triangular factor L, row-major (upper triangle zero).
Definition at line 1902 of file decomposition.hpp.
|
private |
The matrix order.
Definition at line 1900 of file decomposition.hpp.
1.9.8