cheatah
Source BigBrain
Source

stdlib/linalg/routines.hpp

1// Copyright (c) 2026 BigBrain LLC. MIT-licensed (see LICENSE).
2// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.
3#pragma once
5/**
6 * @file routines.hpp
7 * @brief cheatah `linalg` — numpy's linear-algebra API on ndarray, with
8 * SIMD-friendly contiguous kernels (a .purr program writes `linalg.solve(A, b)`).
9 *
10 * `import linalg` to use it (auto-links ndarray). Unit tests:
11 * stdlib/tests/linalg_routines_test.cpp; SIMD reporting tested in
12 * stdlib/tests/linalg_smoke_test.cpp. The suite runs under AddressSanitizer (the `asan`
13 * preset) and Valgrind (security/run-valgrind.sh) on every QA-gate run.
14 *
15 * Routines operate on `ndarray::NDArray` (2-D = matrix, 1-D = vector). They mirror
16 * https://numpy.org/doc/stable/reference/routines.linalg.html and are implemented
17 * in routines.cpp (LU w/ partial pivoting, Cholesky, Householder QR, one-sided
18 * Golub–Reinsch SVD, Householder-tridiagonal + QL symmetric eigen, Hessenberg + shifted-QR general
19 * eigen) at -O3 -march=native so the hot loops auto-vectorize.
20 *
21 * SIMD here is pure compiler auto-vectorization (no intrinsics). On a scalar build
22 * (no vector ISA) every routine still returns identical results, just slower — see
23 * simd.hpp's file comment for the full SIMD model, the no-SIMD behavior, and the
24 * compile-time-dispatch limitation.
25 *
26 * @note `n` below is the matrix dimension. The general eigensolvers `eig`/`eigvals`
27 * return a **complex** spectrum (@ref CNDArray) — a real matrix can have
28 * complex conjugate eigenvalue pairs — while the Hermitian solvers
29 * `eigh`/`eigvalsh` return a guaranteed-real spectrum. The LU/SVD-based scalar
30 * routines (`det`/`slogdet`/`cond`/`matrix_rank`) allocate scratch O(n²) for the
31 * factorization even though they return a scalar; the products and reductions
32 * (`dot`/`matmul`/`trace`/`norm`/…) read their operands in place — zero-copy
33 * when contiguous, packing a strided view once.
34 *
35 * @concurrency every routine is DELIBERATELY single-threaded (no hidden thread pool —
36 * the fastest-per-core contract); the caller composes parallelism across
37 * independent problems.
38 */
39#include <complex>
40#include <vector>
42#include "backend.hpp"
43#include "concepts.hpp"
44#include "enums.hpp"
45#include "ndarray.hpp"
47namespace cheatah::linalg {
49// The routines operate on cheatah::ndarray::NDArray, re-exported unqualified for brevity in the
50// signatures below. The directive below hides this re-export from the API doc generator so it
51// does not emit a phantom duplicate cheatah::linalg::NDArray class in the namespace/XML structure.
52/// @cond INTERNAL
53using ndarray::NDArray;
54/// @endcond
56/// A complex scalar (`std::complex<double>`) — the element type of @ref CNDArray and
57/// the return type of the complex inner products @ref dot / @ref vdot.
58using Cplx = std::complex<double>;
60/// A complex array (`basic_ndarray<std::complex<double>>`) — what the general
61/// eigensolvers return, since a real matrix can have complex eigenvalues. Prints
62/// element-wise as `a+bj` via @ref cheatah::ndarray::to_string.
63using CNDArray = ndarray::basic_ndarray<Cplx>;
65// ---- Matrix and vector products ----
66// Dot / vdot / inner — the scalar reductions — live in backend.hpp as the scalar-out kernel
67// pattern: an allocating front `T dot(a, b)` plus a `void dot(out, a, b)` kernel split by the
68// HostArray/DeviceArray concepts (one pair serving real, complex, host, and — via a device
69// extension — device operands). Same for `trace`.
70// Outer — both the allocating front `outer(a,b)` and the out-parameter kernel `outer(out,a,b)`
71// are the two-layer overload pair in backend.hpp (the matmul pattern).
72// Matmul — both the allocating front `matmul(a,b)` and the out-parameter kernel `matmul(out,a,b)`
73// are the two-layer `template <Field T, template<typename> class Array>` overloads in backend.hpp
74// (one pair serving real, complex, host, and — via a device extension — device operands).
76// ---- complex products (complex inner-product spaces) ----
77// dot / vdot / inner for complex operands are the SAME two-layer templates in backend.hpp,
78// instantiated at T = std::complex<double>; vdot's conjugation is an `if constexpr` branch. No
79// separate symbols. Complex matmul and conj_transpose are likewise the one generic template each
80// in backend.hpp, instantiating the complex element type — no separate complex symbols.
82/// @cond INTERNAL — the HOST kernels of the factorization/solver seam. Every routine below is the
83/// backend.hpp pattern: an allocating front (inline in this header, `requires NumericArray`) that
84/// validates metadata, allocates via `Array<T>::uninitialized`, and calls the same-named out-param
85/// (or scalar-out) kernel UNQUALIFIED — so these HostArray kernels (defined + instantiated in
86/// routines.cpp) serve host arrays, and a device extension's `requires DeviceArray` overloads are
87/// found by ADL. Declared here, before the fronts, so the fronts' unqualified calls see them.
88/// Each kernel's full documentation sits with its declaration further down this header.
89template <ndarray::Field T, template <typename> class Array>
90 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
91void matrix_power(Array<T>& out, const Array<T>& a, long long n);
92template <ndarray::Field T, template <typename> class Array>
93 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
94void cholesky(Array<T>& out, const Array<T>& a);
95template <ndarray::Field T, template <typename> class Array>
96 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
97void qr(Array<T>& q, Array<T>& r, const Array<T>& a);
98template <ndarray::Field T, template <typename> class Array>
99 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
100void svd(Array<T>& u, Array<T>& s, Array<T>& vh, const Array<T>& a);
101template <ndarray::Field T, template <typename> class Array>
102 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
103void svdvals(Array<T>& out, const Array<T>& a);
104template <ndarray::Field T, template <typename> class Array>
105 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
106void eig(Array<ndarray::complex_of_t<T>>& values, Array<ndarray::complex_of_t<T>>& vectors,
107 const Array<T>& a);
108template <ndarray::Field T, template <typename> class Array>
109 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
110void eigvals(Array<ndarray::complex_of_t<T>>& out, const Array<T>& a);
111template <ndarray::Field T, template <typename> class Array>
112 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
113void eigh(Array<ndarray::real_base_t<T>>& values, Array<T>& vectors, const Array<T>& a);
114template <ndarray::Field T, template <typename> class Array>
115 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
116void eigvalsh(Array<ndarray::real_base_t<T>>& out, const Array<T>& a);
117template <ndarray::Field T, template <typename> class Array>
118 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
119void solve(Array<T>& out, const Array<T>& a, const Array<T>& b);
120template <ndarray::Field T, template <typename> class Array>
121 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
122void lstsq(Array<T>& out, const Array<T>& a, const Array<T>& b);
123template <ndarray::Field T, template <typename> class Array>
124 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
125void inv(Array<T>& out, const Array<T>& a);
126template <ndarray::Field T, template <typename> class Array>
127 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
128void pinv(Array<T>& out, const Array<T>& a);
129template <ndarray::Field T, template <typename> class Array>
130 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
131void det(T& out, const Array<T>& a);
132template <ndarray::Field T, template <typename> class Array>
133 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
134void cond(T& out, const Array<T>& a);
135template <ndarray::Field T, template <typename> class Array>
136 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
137void matrix_rank(long long& out, const Array<T>& a);
138template <ndarray::Field T, template <typename> class Array>
139 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
140void norm(T& out, const Array<T>& a);
141/// @endcond
142/**
143 * Integer matrix power Aⁿ (negative n via @ref inv).
144 *
145 * Requires a square matrix (throws otherwise); n == 0 returns the identity, and
146 * negative n first inverts @p a via @ref inv (so it inherits @ref inv's
147 * singular-matrix behavior) before raising to |n|.
148 * @param a square matrix.
149 * @param n exponent.
150 * @return Aⁿ.
151 * @complexity O(n³·log|n|) by binary exponentiation.
152 * @alloc allocates a new NDArray result; the binary-exponentiation @ref matmul steps
153 * allocate their own intermediates.
154 * @test LinalgRoutines.MatrixPower
155 * @crtest LinalgCompileRun.MatrixPower
156 * @systest StdlibE2E.Linalg
157 */
158template <ndarray::Field T, template <typename> class Array>
159 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
160[[nodiscard]] Array<T> matrix_power(const Array<T>& a, long long n) {
161 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
162 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
163 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
164 matrix_power(out, a, n);
165 return out;
167/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
168/**
169 * Matrix power into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
170 * @ref matrix_power (the HOST kernel of the seam pattern; a device extension supplies its own
171 * `requires DeviceArray` overload, found by ADL).
172 * @param out destination; a contiguous n×n matrix, overwritten with Aⁿ.
173 * @param a square matrix.
174 * @param n exponent.
175 * @complexity O(n³·log|n|).
176 * @alloc reuses @p out; the binary-exponentiation products allocate their own scratch.
177 * @test LinalgRoutines.FactorizationOutReusesBuffer
178 */
179template <ndarray::Field T, template <typename> class Array>
180 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
181void matrix_power(Array<T>& out, const Array<T>& a, long long n);
182/// @endcond
183/**
184 * Kronecker product.
185 *
186 * Requires both operands to be 2-D (throws otherwise) and replaces each entry of
187 * @p a with that scalar times the whole of @p b, giving the (m·p)×(k·q) block
188 * matrix; no dimension matching is needed.
189 * @param a m×k matrix.
190 * @param b p×q matrix.
191 * @return (m·p)×(k·q) block product.
192 * @complexity O(n⁴) in the output area.
193 * @alloc allocates a new NDArray result; a non-contiguous operand is packed once into scratch.
194 * @test LinalgRoutines.VdotInnerOuterKron
195 * @crtest LinalgCompileRun.Kron
196 * @systest StdlibE2E.Linalg
197 */
198template <ndarray::Field T, template <typename> class Array>
199 requires NumericArray<Array<T>>
200[[nodiscard]] Array<T> kron(const Array<T>& a, const Array<T>& b);
201/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
202/**
203 * Kronecker product into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
204 * @ref kron, writing the block product straight into @p out with no allocation.
205 * @param out destination; a contiguous (m·p)×(k·q) matrix, overwritten. Must NOT alias @p a or @p b.
206 * @param a m×k matrix.
207 * @param b p×q matrix.
208 * @complexity O(n⁴) in the output area.
209 * @alloc none for contiguous operands (block product written straight into @p out); a
210 * non-contiguous operand is packed once into scratch.
211 * @test LinalgRoutines.KronIntoReusesBuffer
212 */
213template <ndarray::Field T, template <typename> class Array>
214 requires HostArray<Array<T>>
215void kron(Array<T>& out, const Array<T>& a, const Array<T>& b);
216/// @endcond
218// ---- Decompositions ----
219/**
220 * Cholesky factor of a symmetric positive-definite matrix (throws otherwise).
221 *
222 * Requires a square matrix and computes L column by column reading only the
223 * lower triangle of @p a; if any pivot (the diagonal under the square root) is
224 * non-positive it throws "matrix is not positive-definite", which also catches
225 * non-SPD or non-symmetric input.
226 * @param a square SPD matrix.
227 * @return lower-triangular L with A = L·Lᵀ.
228 * @complexity O(n³).
229 * @alloc allocates a new NDArray result; the factor is computed into O(n²) private
230 * scratch, then copied in.
231 * @test LinalgRoutines.CholeskyAndQR
232 * @crtest LinalgCompileRun.Cholesky
233 * @systest StdlibE2E.Linalg
234 */
235template <ndarray::Field T, template <typename> class Array>
236 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
237[[nodiscard]] Array<T> cholesky(const Array<T>& a) { // lower-triangular L (A = L Lᵀ)
238 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
239 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
240 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
241 cholesky(out, a);
242 return out;
244/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
245/**
246 * Cholesky factor into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
247 * @ref cholesky (the HOST kernel of the seam pattern; a device extension supplies its own
248 * `requires DeviceArray` overload, found by ADL).
249 * @param out destination; a contiguous n×n matrix, overwritten with the lower-triangular L.
250 * @param a square SPD matrix.
251 * @complexity O(n³).
252 * @alloc reuses @p out (the factor is computed into private scratch, then copied in).
253 * @test LinalgRoutines.FactorizationOutReusesBuffer
254 */
255template <ndarray::Field T, template <typename> class Array>
256 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
257void cholesky(Array<T>& out, const Array<T>& a);
258/// @endcond
259/** Result of qr(): A = q·r with orthonormal q and upper-triangular r. */
260/// The factors of a QR decomposition. Templated over the container type so a host `qr` yields
261/// `QR<NDArray>` and a device `qr` yields `QR<device_array<T>>`; `ArrT` defaults to `NDArray`
262/// (the host double result), so plain `QR` still names the common host type.
263template <class ArrT = NDArray>
264struct QR {
265 ArrT q; ///< Orthonormal columns, m×n (the Q in A = Q·R).
266 ArrT r; ///< Upper-triangular factor, n×n (the R in A = Q·R).
267};
268/**
269 * Reduced QR via Householder reflections (requires rows ≥ cols).
270 *
271 * Applies successive Householder reflectors to triangularize @p a, returning the
272 * thin/reduced factors; throws "qr requires rows >= cols" for wide matrices.
273 * Rank-deficient columns (zero pivot norm) are skipped, leaving the
274 * corresponding R entries zero.
275 * @param a m×n matrix.
276 * @return @ref QR with q (m×n, orthonormal cols) and r (n×n, upper-triangular).
277 * @complexity O(n³).
278 * @alloc allocates both members; the factorization works in O(m²) private scratch
279 * (a full m×m Q workspace), then copies the reduced factors in.
280 * @test LinalgRoutines.CholeskyAndQR
281 * @crtest LinalgCompileRun.Qr
282 * @systest StdlibE2E.Linalg
283 */
284template <ndarray::Field T, template <typename> class Array>
285 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
286[[nodiscard]] QR<Array<T>> qr(const Array<T>& a) {
287 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
288 const std::size_t m = a.shape()[0], n = a.shape()[1];
289 if (m < n) throw std::runtime_error("linalg: qr requires rows >= cols");
290 QR<Array<T>> out{Array<T>::uninitialized({m, n}), Array<T>::uninitialized({n, n})};
291 qr(out.q, out.r, a);
292 return out;
294/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
295/**
296 * Reduced QR into the caller's buffers (outs FIRST) — the buffer-reuse overload of @ref qr,
297 * filling @p q and @p r instead of allocating a @ref QR (the HOST kernel of the seam pattern; a
298 * device extension supplies its own `requires DeviceArray` overload, found by ADL).
299 * @param q destination for the orthonormal factor; a contiguous m×n matrix, overwritten.
300 * @param r destination for the upper-triangular factor; a contiguous n×n matrix, overwritten.
301 * @param a m×n matrix.
302 * @complexity O(n³).
303 * @alloc reuses @p q and @p r (the factors are computed into private scratch, then copied in).
304 * @test LinalgRoutines.DecompositionOutReusesBuffer
305 */
306template <ndarray::Field T, template <typename> class Array>
307 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
308void qr(Array<T>& q, Array<T>& r, const Array<T>& a);
309/// @endcond
310/** Result of svd(): A = u·diag(s)·vh. */
311/// The factors of a singular value decomposition, templated over the container type (`ArrT`
312/// defaults to `NDArray`, the host double result, so plain `SVD` names the common host type).
313template <class ArrT = NDArray>
314struct SVD {
315 ArrT u; ///< Left singular vectors, m×n.
316 ArrT s; ///< Singular values in descending order (length n).
317 ArrT vh; ///< Right singular vectors transposed, n×n (the Vᵀ in A = u·diag(s)·Vᵀ).
318};
319/**
320 * Singular value decomposition (Golub–Reinsch; requires rows ≥ cols).
321 *
322 * Reduces @p a to upper-bidiagonal form by Householder reflections, then diagonalizes
323 * it with implicit-shift QR (accumulating U and V), and sorts the singular values
324 * descending — the world-standard dense SVD (what LAPACK's dgesvd reduces to). Throws
325 * "svd requires rows >= cols" for wide matrices (transpose first); singular values come
326 * out non-negative.
327 * @param a m×n matrix.
328 * @return @ref SVD with u (m×n), s (descending singular values), vh (n×n = Vᵀ).
329 * @complexity iterative O(n³).
330 * @alloc allocates all members; the Golub–Reinsch reduction allocates its own O(m·n) workspace.
331 * @test LinalgRoutines.SvdAndEigh
332 * @crtest LinalgCompileRun.Svd
333 * @systest StdlibE2E.Linalg
334 */
335template <ndarray::Field T, template <typename> class Array>
336 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
337[[nodiscard]] SVD<Array<T>> svd(const Array<T>& a) {
338 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
339 const std::size_t m = a.shape()[0], n = a.shape()[1];
340 if (m < n) throw std::runtime_error("linalg: svd requires rows >= cols (transpose otherwise)");
341 SVD<Array<T>> out{Array<T>::uninitialized({m, n}), Array<T>::uninitialized({n}),
342 Array<T>::uninitialized({n, n})};
343 svd(out.u, out.s, out.vh, a);
344 return out;
346/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
347/**
348 * Full SVD into the caller's buffers (outs FIRST) — the buffer-reuse overload of @ref svd,
349 * filling @p u, @p s and @p vh instead of allocating an @ref SVD (the HOST kernel of the seam
350 * pattern; a device extension supplies its own `requires DeviceArray` overload, found by ADL).
351 * @param u destination for the left singular vectors; a contiguous m×n matrix, overwritten.
352 * @param s destination for the singular values; a contiguous length-n vector, overwritten.
353 * @param vh destination for Vᵀ; a contiguous n×n matrix, overwritten.
354 * @param a m×n matrix (rows ≥ cols).
355 * @complexity iterative O(n³).
356 * @alloc reuses @p u, @p s, @p vh (the factors are computed into private scratch, then copied in).
357 * @test LinalgRoutines.DecompositionOutReusesBuffer
358 */
359template <ndarray::Field T, template <typename> class Array>
360 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
361void svd(Array<T>& u, Array<T>& s, Array<T>& vh, const Array<T>& a);
362/// @endcond
363/**
364 * Singular values only (≈ `numpy.linalg.svd(a, compute_uv=False)` / `svdvals`).
365 *
366 * Runs the same Golub–Reinsch reduction as @ref svd but takes the **values-only** fast
367 * path — it never accumulates U or V, and skips the (dominant) U/V Givens rotations in
368 * the QR sweep — so it is several times faster than the full decomposition. Accepts any
369 * shape (singular values of `a` and `aᵀ` coincide).
370 * @param a m×n matrix.
371 * @return length-min(m,n) vector of singular values, descending.
372 * @complexity iterative O(n³), but a large constant factor below @ref svd.
373 * @alloc allocates a new NDArray result; the reduction still allocates its O(m·n)
374 * workspace (the values-only path skips the U/V accumulation work and result
375 * copies, not the working buffers).
376 * @test LinalgRoutines.SvdAndEigh
377 * @crtest LinalgCompileRun.Svdvals
378 * @systest StdlibE2E.Linalg
379 */
380template <ndarray::Field T, template <typename> class Array>
381 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
382[[nodiscard]] Array<T> svdvals(const Array<T>& a) {
383 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
384 const std::size_t m = a.shape()[0], n = a.shape()[1];
385 Array<T> out = Array<T>::uninitialized({m < n ? m : n});
386 svdvals(out, a);
387 return out;
389/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
390/**
391 * Singular values into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
392 * @ref svdvals (the HOST kernel of the seam pattern; a device extension supplies its own
393 * `requires DeviceArray` overload, found by ADL).
394 * @param out destination; a contiguous length-min(m,n) vector, overwritten with the descending values.
395 * @param a m×n matrix.
396 * @complexity iterative O(n³).
397 * @alloc reuses @p out; the Golub–Reinsch reduction allocates its own scratch.
398 * @test LinalgRoutines.FactorizationOutReusesBuffer
399 */
400template <ndarray::Field T, template <typename> class Array>
401 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
402void svdvals(Array<T>& out, const Array<T>& a);
403/// @endcond
405// ---- Matrix eigenvalues ----
406/** Result of eigh(): a real spectrum — column j of vectors is the eigenvector for values[j]. */
407template <class ArrT = NDArray>
408struct Eig {
409 ArrT values; ///< Eigenvalues (length n), real.
410 ArrT vectors; ///< Eigenvectors as columns: column j matches values[j] (empty if not computed).
411};
412/**
413 * Result of the general eig(): a **complex** spectrum, since a real matrix can have
414 * complex conjugate eigenvalue pairs. Column j of vectors is the eigenvector for values[j].
415 */
416template <class ArrT = CNDArray>
417struct EigC {
418 ArrT values; ///< Eigenvalues (length n), complex.
419 ArrT vectors; ///< Eigenvectors as columns: column j matches values[j].
420};
421/// Result of the complex Hermitian eigh(): **real** eigenvalues with **complex** eigenvectors —
422/// so the two members have DIFFERENT container/element types (`ValsT` real, `VecsT` complex).
423/// Defaults `NDArray`/`CNDArray` are the host result, so plain `EighC` names the common host type.
424template <class ValsT = NDArray, class VecsT = CNDArray>
425struct EighC {
426 ValsT values; ///< Eigenvalues (length n), real and descending.
427 VecsT vectors; ///< Eigenvectors as columns: column j matches values[j].
428};
429/// The result type of the unified @ref eigh: for a real element, `Eig<Array<T>>` (real values +
430/// vectors); for a complex element, `EighC<Array<real>, Array<T>>` (real values, complex vectors).
431/// eigh's return type differs per element, so it is expressed here (not `auto`) so the header
432/// declaration knows it without seeing the definition.
433template <ndarray::Field T, template <typename> class Array>
434using eigh_result_t = std::conditional_t<ndarray::is_complex_v<T>,
435 EighC<Array<ndarray::real_base_t<T>>, Array<T>>,
436 Eig<Array<T>>>;
437/**
438 * Eigen-decomposition of a general square matrix (**complex** spectrum and
439 * eigenvectors).
440 *
441 * For a symmetric @p a it delegates to @ref eigh (promoted to complex with zero
442 * imaginary part); otherwise it uses Hessenberg reduction + shifted QR for the
443 * eigenvalues, then **inverse iteration** for each eigenvector. A real matrix with a
444 * complex conjugate pair (e.g. a rotation) yields those complex eigenvalues and
445 * eigenvectors rather than throwing. Throws on a non-square matrix or if the QR
446 * iteration fails to converge.
447 * @param a square matrix.
448 * @return @ref EigC with complex values and matching complex eigenvector columns.
449 * @complexity iterative O(n³) via Hessenberg + shifted QR for the eigenvalues; the
450 * general (non-symmetric) eigenvectors add O(n⁴) — one inverse iteration, each
451 * with its own O(n³) complex LU factorization, per eigenvalue (a symmetric @p a
452 * stays O(n³) via @ref eigh, which accumulates the vectors in the QL sweep).
453 * @alloc allocates both members; plus O(n²) factorization scratch (on the non-symmetric
454 * path, a fresh complex n×n LU per eigenvalue).
455 * @test LinalgRoutines.GeneralEig
456 * @crtest LinalgCompileRun.Eig
457 * @systest StdlibE2E.Linalg
458 */
459template <ndarray::Field T, template <typename> class Array>
460 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
461[[nodiscard]] EigC<Array<ndarray::complex_of_t<T>>> eig(const Array<T>& a) { // general square matrix
462 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
463 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
464 const std::size_t n = a.shape()[0];
465 using C = Array<ndarray::complex_of_t<T>>;
466 EigC<C> out{C::uninitialized({n}), C::uninitialized({n, n})};
467 eig(out.values, out.vectors, a);
468 return out;
470/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
471/**
472 * General eigendecomposition into the caller's buffers (outs FIRST) — the buffer-reuse overload of
473 * @ref eig, filling @p values and @p vectors instead of allocating an @ref EigC (the HOST kernel
474 * of the seam pattern; a device extension supplies its own `requires DeviceArray` overload).
475 * @param values destination for the complex eigenvalues; a contiguous length-n vector, overwritten.
476 * @param vectors destination for the complex eigenvectors (columns); a contiguous n×n matrix, overwritten.
477 * @param a square matrix.
478 * @complexity iterative O(n³) for the eigenvalues; general (non-symmetric) eigenvectors
479 * add O(n⁴) (inverse iteration per eigenvalue — see @ref eig).
480 * @alloc reuses @p values and @p vectors (computed into private scratch, then copied in).
481 * @test LinalgRoutines.DecompositionOutReusesBuffer
482 */
483template <ndarray::Field T, template <typename> class Array>
484 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
485void eig(Array<ndarray::complex_of_t<T>>& values, Array<ndarray::complex_of_t<T>>& vectors,
486 const Array<T>& a);
487/// @endcond
488/**
489 * Eigenvalues of a general square matrix (**complex**), descending.
490 *
491 * Routes symmetric input through tridiagonal QL and everything else through
492 * Hessenberg + shifted QR, then sorts the result descending (by real part, then by
493 * imaginary part). A real matrix with a complex conjugate pair yields those complex
494 * eigenvalues rather than throwing. Throws on a non-square matrix or non-convergence
495 * of the QR iteration.
496 * @param a square matrix.
497 * @return length-n complex vector of eigenvalues.
498 * @complexity iterative O(n³).
499 * @alloc allocates a new CNDArray result; the iteration allocates its own O(n²) scratch.
500 * @test LinalgRoutines.SvdAndEigh
501 * @crtest LinalgCompileRun.Eigvals
502 * @systest StdlibE2E.Linalg
503 */
504template <ndarray::Field T, template <typename> class Array>
505 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
506[[nodiscard]] Array<ndarray::complex_of_t<T>> eigvals(const Array<T>& a) {
507 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
508 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
509 Array<ndarray::complex_of_t<T>> out =
510 Array<ndarray::complex_of_t<T>>::uninitialized({a.shape()[0]});
511 eigvals(out, a);
512 return out;
514/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
515/**
516 * General eigenvalues into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
517 * @ref eigvals (the HOST kernel of the seam pattern; a device extension supplies its own
518 * `requires DeviceArray` overload, found by ADL).
519 * @param out destination; a contiguous length-n complex vector, overwritten with the descending spectrum.
520 * @param a square matrix.
521 * @complexity iterative O(n³).
522 * @alloc reuses @p out; the Hessenberg + shifted-QR iteration allocates its own scratch.
523 * @test LinalgRoutines.FactorizationOutReusesBuffer
524 */
525template <ndarray::Field T, template <typename> class Array>
526 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
527void eigvals(Array<ndarray::complex_of_t<T>>& out, const Array<T>& a);
528/// @endcond
529/**
530 * Eigen-decomposition of a symmetric matrix (Householder tridiagonalization + QL).
531 *
532 * Reduces @p a to tridiagonal form by Householder reflections, then diagonalizes it
533 * with implicit-shift QL, returning real eigenvalues sorted descending with matching
534 * eigenvector columns; it reads the full matrix and assumes symmetry rather than
535 * checking it, so asymmetric input yields meaningless results. Throws on a non-square
536 * matrix (or if the QL iteration fails to converge).
537 * @param a square symmetric matrix.
538 * @return @ref Eig with descending values and matching eigenvectors.
539 * @complexity iterative O(n³).
540 * @alloc allocates both members; the solver allocates its own O(n²) scratch (a complex
541 * Hermitian input first embeds into a 2n×2n real matrix).
542 * @test LinalgRoutines.SvdAndEigh
543 * @test LinalgRoutines.ComplexHermitianEigh
544 * @crtest LinalgCompileRun.Eigh
545 * @crtest LinalgCompileRun.EighComplex
546 * @systest StdlibE2E.Linalg
547 * @systest StdlibE2E.LinalgComplex
548 */
549template <ndarray::Field T, template <typename> class Array>
550 requires NumericArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
551[[nodiscard]] eigh_result_t<T, Array> eigh(const Array<T>& a) { // symmetric / Hermitian
552 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
553 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
554 const std::size_t n = a.shape()[0];
555 eigh_result_t<T, Array> out{Array<ndarray::real_base_t<T>>::uninitialized({n}),
556 Array<T>::uninitialized({n, n})};
557 eigh(out.values, out.vectors, a);
558 return out;
560/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
561/**
562 * Symmetric/Hermitian eigendecomposition into the caller's buffers (outs FIRST) — the buffer-reuse
563 * overload of @ref eigh, filling @p values and @p vectors instead of allocating a result struct
564 * (the HOST kernel of the seam pattern; a device extension supplies its own `requires DeviceArray`
565 * overload). ONE two-layer kernel serving the real symmetric AND complex Hermitian paths: values
566 * are always the real spectrum, vectors match the input element.
567 * @param values destination for the real eigenvalues; a contiguous length-n vector, overwritten.
568 * @param vectors destination for the eigenvectors (columns); a contiguous n×n matrix, overwritten.
569 * @param a square symmetric (real) / Hermitian (complex) matrix.
570 * @complexity iterative O(n³).
571 * @alloc reuses @p values and @p vectors (computed into private scratch, then copied in).
572 * @test LinalgRoutines.DecompositionOutReusesBuffer
573 */
574template <ndarray::Field T, template <typename> class Array>
575 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
576void eigh(Array<ndarray::real_base_t<T>>& values, Array<T>& vectors, const Array<T>& a);
577/// @endcond
578/**
579 * Eigenvalues of a symmetric matrix, descending (tridiagonal QL).
580 *
581 * Same tridiagonalization + QL as @ref eigh but **skips the eigenvector accumulation
582 * entirely** (the bulk of the work), so it is roughly twice as fast as `eigh`; assumes
583 * (does not verify) symmetry and throws on a non-square matrix.
584 * ONE two-layer template collapsing the former real and complex (Hermitian) overloads: a real
585 * element takes the symmetric path, a complex element the Hermitian path (`if constexpr`). The
586 * spectrum is always REAL, returned as `Array<real_base_t<T>>`.
587 * @tparam T the element type (`double` or `std::complex<double>`); @tparam Array the container.
588 * @param a square symmetric (real) / Hermitian (complex) matrix.
589 * @return length-n vector of real eigenvalues.
590 * @complexity iterative O(n³).
591 * @alloc allocates a new result; the solver allocates its own O(n²) scratch (a complex
592 * Hermitian input first embeds into a 2n×2n real matrix).
593 * @test LinalgRoutines.EigvalshSymmetric
594 * @test LinalgRoutines.ComplexHermitianEigh
595 * @crtest LinalgCompileRun.Eigvalsh
596 * @crtest LinalgCompileRun.EigvalshComplex
597 * @systest StdlibE2E.Linalg
598 * @systest StdlibE2E.LinalgComplex
599 */
600template <ndarray::Field T, template <typename> class Array>
601 requires NumericArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
602[[nodiscard]] Array<ndarray::real_base_t<T>> eigvalsh(const Array<T>& a) {
603 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
604 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
605 Array<ndarray::real_base_t<T>> out =
606 Array<ndarray::real_base_t<T>>::uninitialized({a.shape()[0]});
607 eigvalsh(out, a);
608 return out;
610/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
611/**
612 * Symmetric/Hermitian eigenvalues into the caller's buffer @p out (out FIRST) — the buffer-reuse
613 * overload of @ref eigvalsh (the HOST kernel of the seam pattern; a device extension supplies its
614 * own `requires DeviceArray` overload). One two-layer kernel: the complex Hermitian path is the
615 * same template at T = std::complex<double>, still writing the REAL spectrum.
616 * @param out destination; a contiguous length-n vector, overwritten with the descending eigenvalues.
617 * @param a square symmetric (real) / Hermitian (complex) matrix.
618 * @complexity iterative O(n³).
619 * @alloc reuses @p out; the tridiagonal-QL solver allocates its own scratch.
620 * @test LinalgRoutines.FactorizationOutReusesBuffer
621 */
622template <ndarray::Field T, template <typename> class Array>
623 requires HostArray<Array<T>> && ndarray::FloatingPoint<ndarray::real_base_t<T>>
624void eigvalsh(Array<ndarray::real_base_t<T>>& out, const Array<T>& a);
625/// @endcond
626// EighC (real values + complex vectors) and the unified two-layer `eigh` are declared above with
627// the other eig-family structs; the complex Hermitian eigh is that template at T = complex<double>,
628// returning EighC<NDArray, CNDArray> via the if-constexpr Hermitian branch — and its buffer-reuse
629// form is the SAME two-layer eigh kernel above at T = complex<double> (values NDArray&, vectors
630// CNDArray&). Likewise complex eigvalsh is the eigvalsh kernel at T = complex<double>.
631// ---- Norms and other numbers ----
632/**
633 * Norm: L2 for vectors, Frobenius for matrices.
634 *
635 * Dispatches on rank: 1-D (or lower) inputs get the Euclidean L2 norm, 2-D
636 * inputs the Frobenius norm; either way it is the square root of the sum of
637 * squared entries. Two-layer over the element and container like every routine (the scalar-out
638 * kernel `norm(out, a)` is the host/device seam; host `double` is the shipped instantiation).
639 * @param a vector or matrix.
640 * @return √Σ xᵢ².
641 * @complexity O(n) for vectors / O(n²) for matrices.
642 * @alloc none for a contiguous operand (summed in place); a non-contiguous view packs
643 * once. Returns a double.
644 * @test LinalgRoutines.NormAndRank
645 * @crtest LinalgCompileRun.Norm
646 * @systest StdlibE2E.Linalg
647 */
648template <ndarray::Field T, template <typename> class Array>
649 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
650[[nodiscard]] T norm(const Array<T>& a) { // default: Frobenius / L2
651 T out;
652 norm(out, a);
653 return out;
655/**
656 * 2-norm condition number σ_max/σ_min (∞ if singular).
657 *
658 * Takes the ratio of largest to smallest singular value from a Golub–Reinsch SVD
659 * (transposing internally for wide matrices); returns +infinity when the
660 * smallest singular value is exactly zero (singular/rank-deficient).
661 * @param a matrix.
662 * @return condition number.
663 * @complexity iterative O(n³) via SVD.
664 * @alloc allocates scratch O(n²) for the factorization; returns a double.
665 * @test LinalgRoutines.SlogdetAndCond
666 * @crtest LinalgCompileRun.Cond
667 * @systest StdlibE2E.Linalg
668 */
669template <ndarray::Field T, template <typename> class Array>
670 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
671[[nodiscard]] T cond(const Array<T>& a) {
672 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
673 T out;
674 cond(out, a);
675 return out;
677/**
678 * Determinant via LU with partial pivoting.
679 *
680 * Computes the product of the LU pivots times the permutation sign; requires a
681 * square matrix (throws otherwise). A singular matrix yields a determinant of
682 * (or extremely near) zero rather than an error.
683 * @param a square matrix.
684 * @return det(A).
685 * @complexity O(n³).
686 * @alloc allocates scratch O(n²) for the factorization; returns a double.
687 * @test LinalgRoutines.SolveDetInv
688 * @crtest LinalgCompileRun.Det
689 * @systest StdlibE2E.Linalg
690 */
691template <ndarray::Field T, template <typename> class Array>
692 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
693[[nodiscard]] T det(const Array<T>& a) {
694 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
695 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
696 T out;
697 det(out, a);
698 return out;
700/**
701 * Numerical rank from SVD singular-value thresholding.
702 *
703 * Counts singular values above a tolerance scaled by the largest singular value
704 * and the matrix size (the standard numpy-style threshold); accepts any shape,
705 * transposing wide matrices internally.
706 * @param a matrix.
707 * @return rank.
708 * @complexity iterative O(n³) via SVD.
709 * @alloc allocates scratch O(n²) for the factorization.
710 * @test LinalgRoutines.NormAndRank
711 * @crtest LinalgCompileRun.MatrixRank
712 * @systest StdlibE2E.Linalg
713 */
714template <ndarray::Field T, template <typename> class Array>
715 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
716[[nodiscard]] long long matrix_rank(const Array<T>& a) {
717 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
718 long long out;
719 matrix_rank(out, a);
720 return out;
722/** Result of slogdet(): det(A) = sign·exp(logabsdet). */
723struct SLogDet {
724 double sign; ///< Sign of the determinant (−1, 0, or +1).
725 double logabsdet; ///< Natural log of |det(A)|, so det(A) = sign·exp(logabsdet).
726};
727/// @cond INTERNAL — the scalar-out HOST kernel of the seam pattern (declared after @ref SLogDet,
728/// which its out-parameter needs; a device extension supplies its own DeviceArray overload).
729template <ndarray::Field T, template <typename> class Array>
730 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
731void slogdet(SLogDet& out, const Array<T>& a);
732/// @endcond
733/**
734 * Sign and log|det| via LU (overflow-safe determinant).
735 *
736 * Sums the logs of the absolute LU pivots (avoiding the over/underflow of a raw
737 * product) and tracks the sign from the pivot signs and permutation parity;
738 * requires a square matrix (throws otherwise). A singular matrix gives a hugely
739 * negative logabsdet rather than −infinity, since a zero pivot is nudged to a
740 * tiny value during factorization.
741 * @param a square matrix.
742 * @return @ref SLogDet.
743 * @complexity O(n³).
744 * @alloc allocates scratch O(n²) for the factorization (the struct members are plain doubles).
745 * @test LinalgRoutines.SlogdetAndCond
746 * @crtest LinalgCompileRun.Slogdet
747 * @systest StdlibE2E.Linalg
748 */
749template <ndarray::Field T, template <typename> class Array>
750 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
751[[nodiscard]] SLogDet slogdet(const Array<T>& a) {
752 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
753 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
754 SLogDet out{};
755 slogdet(out, a);
756 return out;
758// Trace — the allocating front `trace(a)` and the scalar-out kernel `trace(out, a)` are the
759// backend.hpp reduction pattern (host kernel here in routines.cpp; a device extension adds its
760// own DeviceArray overload, found by ADL).
762// ---- Solving equations and inverting matrices ----
763/**
764 * Solve A·x = b via LU with partial pivoting.
765 *
766 * Factorizes @p a once then does forward/back substitution against @p b;
767 * requires @p a square and @p b a vector of matching length (throws otherwise).
768 * A singular @p a does not throw but yields a garbage/overflowing solution
769 * (pivots are nudged off zero rather than detected).
770 * @param a square coefficient matrix.
771 * @param b right-hand-side vector.
772 * @return solution x.
773 * @complexity O(n³).
774 * @alloc allocates a new NDArray result; the LU factorization allocates its own O(n²) scratch.
775 * @test LinalgRoutines.SolveDetInv
776 * @crtest LinalgCompileRun.Solve
777 * @systest StdlibE2E.Linalg
778 */
779template <ndarray::Field T, template <typename> class Array>
780 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
781[[nodiscard]] Array<T> solve(const Array<T>& a, const Array<T>& b) { // A x = b
782 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
783 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
784 if (vector_len(b) != a.shape()[0]) throw std::runtime_error("linalg: solve dimension mismatch");
785 Array<T> out = Array<T>::uninitialized({a.shape()[0]});
786 solve(out, a, b);
787 return out;
789/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
790/**
791 * Solve into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of @ref solve.
792 * @param out destination; a contiguous length-n vector, overwritten with the solution x.
793 * @param a square coefficient matrix.
794 * @param b right-hand-side vector.
795 * @complexity O(n³).
796 * @alloc reuses @p out; the LU factorization allocates its own scratch.
797 * @test LinalgRoutines.FactorizationOutReusesBuffer
798 */
799template <ndarray::Field T, template <typename> class Array>
800 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
801void solve(Array<T>& out, const Array<T>& a, const Array<T>& b);
802/// @endcond
803/**
804 * Least-squares solution min‖A·x − b‖ (computed as @ref pinv (a)·b).
805 *
806 * Forms the Moore–Penrose pseudo-inverse via SVD and multiplies it by @p b, so
807 * it handles over- and under-determined systems and returns the minimum-norm
808 * solution for rank-deficient @p a; @p b must be conformable for the
809 * @ref matmul step.
810 * @param a m×n matrix.
811 * @param b right-hand side.
812 * @return minimizing x.
813 * @complexity iterative O(n³) via SVD.
814 * @alloc allocates a new NDArray result; plus the intermediate n×m pseudo-inverse and
815 * its SVD scratch.
816 * @test LinalgRoutines.Lstsq
817 * @crtest LinalgCompileRun.Lstsq
818 * @systest StdlibE2E.Linalg
819 */
820template <ndarray::Field T, template <typename> class Array>
821 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
822[[nodiscard]] Array<T> lstsq(const Array<T>& a, const Array<T>& b) { // least-squares solution
823 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
824 if (b.ndim() != 2) throw std::runtime_error("linalg: matmul expects 2-D matrices");
825 if (a.shape()[0] != b.shape()[0])
826 throw std::runtime_error("linalg: matmul inner dimension mismatch");
827 Array<T> out = Array<T>::uninitialized({a.shape()[1], b.shape()[1]});
828 lstsq(out, a, b);
829 return out;
831/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
832/**
833 * Least-squares solution into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
834 * @ref lstsq. Routes through the @ref matmul out-param so the final product is written into @p out
835 * with no allocation.
836 * @param out destination; a contiguous array of the solution's shape, overwritten.
837 * @param a m×n matrix.
838 * @param b right-hand side.
839 * @complexity iterative O(n³).
840 * @alloc reuses @p out; allocates the intermediate n×m pseudo-inverse and its SVD scratch.
841 * @test LinalgRoutines.FactorizationOutReusesBuffer
842 */
843template <ndarray::Field T, template <typename> class Array>
844 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
845void lstsq(Array<T>& out, const Array<T>& a, const Array<T>& b);
846/// @endcond
847/**
848 * Matrix inverse via LU with partial pivoting.
849 *
850 * Factorizes @p a once and back-solves against each identity column; requires a
851 * square matrix (throws otherwise). A singular @p a does not throw but produces
852 * garbage/overflowing entries since zero pivots are nudged rather than detected.
853 * @param a square matrix.
854 * @return A⁻¹.
855 * @complexity O(n³).
856 * @alloc allocates a new NDArray result; the LU factorization and the whole-identity
857 * back-solve allocate their own O(n²) scratch.
858 * @test LinalgRoutines.SolveDetInv
859 * @crtest LinalgCompileRun.Inv
860 * @systest StdlibE2E.Linalg
861 */
862template <ndarray::Field T, template <typename> class Array>
863 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
864[[nodiscard]] Array<T> inv(const Array<T>& a) {
865 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
866 if (a.shape()[0] != a.shape()[1]) throw std::runtime_error("linalg: expected a square matrix");
867 Array<T> out = Array<T>::uninitialized({a.shape()[0], a.shape()[0]});
868 inv(out, a);
869 return out;
871/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
872/**
873 * Inverse into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of @ref inv.
874 * @param out destination; a contiguous n×n matrix, overwritten with A⁻¹.
875 * @param a square matrix.
876 * @complexity O(n³).
877 * @alloc reuses @p out; the LU factorization allocates its own scratch.
878 * @test LinalgRoutines.FactorizationOutReusesBuffer
879 */
880template <ndarray::Field T, template <typename> class Array>
881 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
882void inv(Array<T>& out, const Array<T>& a);
883/// @endcond
884/**
885 * Moore–Penrose pseudo-inverse via SVD (any shape).
886 *
887 * Computes V·diag(1/σ)·Uᵀ from a Golub–Reinsch SVD, transposing wide matrices
888 * internally so any shape works; singular values at or below a size-scaled
889 * tolerance are dropped (treated as zero) so it stays well-defined for
890 * rank-deficient input.
891 * @param a m×n matrix.
892 * @return n×m pseudo-inverse.
893 * @complexity iterative O(n³) via SVD.
894 * @alloc allocates a new NDArray result; the SVD and the assembly allocate their own
895 * O(m·n) scratch.
896 * @test LinalgRoutines.PinvCondRankOnWideMatrix
897 * @crtest LinalgCompileRun.Pinv
898 * @systest StdlibE2E.Linalg
899 */
900template <ndarray::Field T, template <typename> class Array>
901 requires NumericArray<Array<T>> && ndarray::FloatingPoint<T>
902[[nodiscard]] Array<T> pinv(const Array<T>& a) { // Moore–Penrose pseudo-inverse
903 if (a.ndim() != 2) throw std::runtime_error("linalg: expected a 2-D matrix");
904 Array<T> out = Array<T>::uninitialized({a.shape()[1], a.shape()[0]});
905 pinv(out, a);
906 return out;
908/// @cond INTERNAL — the allocation-free out-parameter variant (see README: buffer reuse)
909/**
910 * Pseudo-inverse into the caller's buffer @p out (out FIRST) — the buffer-reuse overload of
911 * @ref pinv.
912 * @param out destination; a contiguous n×m matrix (for an m×n input), overwritten with the pseudo-inverse.
913 * @param a m×n matrix.
914 * @complexity iterative O(n³).
915 * @alloc reuses @p out; the Golub–Reinsch SVD allocates its own scratch.
916 * @test LinalgRoutines.FactorizationOutReusesBuffer
917 */
918template <ndarray::Field T, template <typename> class Array>
919 requires HostArray<Array<T>> && ndarray::FloatingPoint<T>
920void pinv(Array<T>& out, const Array<T>& a);
921/// @endcond
923} // namespace cheatah::linalg