linalg

Both operands are Array<T> (so host⊗device / element mixes fail to deduce and are compile errors); requires both to be 2-D with matching inner dimensions — or both 3-D for the BATCHED product [B,M,K] @ [B,K,N] → [B,M,N] (equal batch counts, strict: no broadcast batching). Allocates the result via Array<T>::uninitialized and fills it through the out-parameter kernel — the host SIMD path, or a device shader when Array is a device container.

Reads only host-resident shape metadata, so it is valid for ANY located container, device arrays included; the shared validation step of every vector front below.

Flattens each operand to a vector (1-D, or 2-D with a size-1 row/column) and throws if either is not vector-shaped or the lengths differ. Both operands are Array<T> (the deduction firewall).

For a REAL element this is the bilinear Σ aᵢbᵢ (identical to dot and inner); for a complex element it is the conjugate-linear Hermitian inner product ⟨a, b⟩ = Σ conj(aᵢ)·bᵢ (numpy's vdot, conjugating the first argument); vdot(a, a) is ‖a‖².

Requires a 2-D matrix (throws otherwise); rectangular matrices sum min(r, c) diagonal entries.

Flattens both operands to vectors and forms the full rank-1 matrix; any pair of vector lengths is accepted (no matching constraint). Allocates the result via Array<T>::uninitialized and fills it through the out-parameter kernel — the host SIMD path, or a device shader when Array is a device container (selected by concept).

A matrix is Hermitian iff conj_transpose(A) == A.

Requires both operands to be 2-D (throws otherwise) and replaces each entry of a with that scalar times the whole of b, giving the (m·p)×(k·q) block matrix; no dimension matching is needed.

Requires a square matrix (throws otherwise); n == 0 returns the identity, and negative n first inverts a via inv (so it inherits inv's singular-matrix behavior) before raising to |n|.

Dispatches on rank: 1-D (or lower) inputs get the Euclidean L2 norm, 2-D inputs the Frobenius norm; either way it is the square root of the sum of squared entries. Two-layer over the element and container like every routine (the scalar-out kernel norm(out, a) is the host/device seam; host double is the shipped instantiation).

Factorizes a once then does forward/back substitution against b; requires a square and b a vector of matching length (throws otherwise). A singular a does not throw but yields a garbage/overflowing solution (pivots are nudged off zero rather than detected).

Computes the product of the LU pivots times the permutation sign; requires a square matrix (throws otherwise). A singular matrix yields a determinant of (or extremely near) zero rather than an error.

Sums the logs of the absolute LU pivots (avoiding the over/underflow of a raw product) and tracks the sign from the pivot signs and permutation parity; requires a square matrix (throws otherwise). A singular matrix gives a hugely negative logabsdet rather than −infinity, since a zero pivot is nudged to a tiny value during factorization.

Factorizes a once and back-solves against each identity column; requires a square matrix (throws otherwise). A singular a does not throw but produces garbage/overflowing entries since zero pivots are nudged rather than detected.

Forms the Moore–Penrose pseudo-inverse via SVD and multiplies it by b, so it handles over- and under-determined systems and returns the minimum-norm solution for rank-deficient a; b must be conformable for the matmul step.

Requires a square matrix and computes L column by column reading only the lower triangle of a; if any pivot (the diagonal under the square root) is non-positive it throws "matrix is not positive-definite", which also catches non-SPD or non-symmetric input.

Applies successive Householder reflectors to triangularize a, returning the thin/reduced factors; throws "qr requires rows >= cols" for wide matrices. Rank-deficient columns (zero pivot norm) are skipped, leaving the corresponding R entries zero.

Reduces a to upper-bidiagonal form by Householder reflections, then diagonalizes it with implicit-shift QR (accumulating U and V), and sorts the singular values descending — the world-standard dense SVD (what LAPACK's dgesvd reduces to). Throws "svd requires rows >= cols" for wide matrices (transpose first); singular values come out non-negative.

Runs the same Golub–Reinsch reduction as svd but takes the values-only fast path — it never accumulates U or V, and skips the U/V Givens rotations in the QR sweep — so it beats the full decomposition by a margin that grows with n (see the linalg benchmark page). Accepts any shape (singular values of a and aᵀ coincide).

Computes V·diag(1/σ)·Uᵀ from a Golub–Reinsch SVD, transposing wide matrices internally so any shape works; singular values at or below a size-scaled tolerance are dropped (treated as zero) so it stays well-defined for rank-deficient input.

Takes the ratio of largest to smallest singular value from a Golub–Reinsch SVD (transposing internally for wide matrices); returns +infinity when the smallest singular value is exactly zero (singular/rank-deficient).

Counts singular values above a tolerance scaled by the largest singular value and the matrix size (the standard numpy-style threshold); accepts any shape, transposing wide matrices internally.

Reduces a to tridiagonal form by Householder reflections, then diagonalizes it with implicit-shift QL, returning real eigenvalues sorted descending with matching eigenvector columns; it reads the full matrix and assumes symmetry rather than checking it, so asymmetric input yields meaningless results. Throws on a non-square matrix (or if the QL iteration fails to converge).

Same tridiagonalization + QL as eigh but skips the eigenvector accumulation entirely, so it runs well under eigh's cost at large n (see the linalg benchmark page); assumes (does not verify) symmetry and throws on a non-square matrix. ONE two-layer template collapsing the former real and complex (Hermitian) overloads: a real element takes the symmetric path, a complex element the Hermitian path (if constexpr). The spectrum is always REAL, returned as Array<real_base_t<T>>.

For a symmetric a it delegates to eigh (promoted to complex with zero imaginary part); otherwise it uses Hessenberg reduction + shifted QR for the eigenvalues, then inverse iteration for each eigenvector. A real matrix with a complex conjugate pair (e.g. a rotation) yields those complex eigenvalues and eigenvectors rather than throwing. Throws on a non-square matrix or if the QR iteration fails to converge.

Routes symmetric input through tridiagonal QL and everything else through Hessenberg + shifted QR, then sorts the result descending (by real part, then by imaginary part). A real matrix with a complex conjugate pair yields those complex eigenvalues rather than throwing. Throws on a non-square matrix or non-convergence of the QR iteration.

"AVX2;FMA" (x86-64), "NEON" (ARM), or "scalar".

Reflects compile-time target flags (e.g. -march=native), not a runtime CPUID probe.

1 = scalar, 2 = SSE2/NEON, 4 = AVX, 8 = AVX-512; useful for sizing blocked kernels.

dot/inner are bilinear (Conj::None, Σ aᵢbᵢ); vdot is the conjugate-linear Hermitian inner product (Conj::Conjugate, Σ conj(aᵢ)·bᵢ). For a REAL element the conjugate is the identity, so both fold to the same code via if constexpr (is_complex_v<T> && mode == Conj::Conjugate).

Undefined for a type with no value_type (which simply makes the concepts below unsatisfied for it, never a hard error).

Prints element-wise as a+bj via cheatah::ndarray::to_string.

eigh's return type differs per element, so it is expressed here (not auto) so the header declaration knows it without seeing the definition.