The NumericAL public API. More...

Namespaces
namespace	blas

namespace	detail

Classes
class	cholesky_decomposition
	A Cholesky factorisation `A = L Lᴴ` of a Hermitian positive-definite matrix. More...

class	convergence_error
	Thrown when an iterative algorithm has not converged within its bound of sweeps. More...

class	dimension_error
	Thrown when an operation's operands have incompatible shapes. More...

class	eigen_decomposition
	The eigenvalues, and optionally the right eigenvectors, of a general square matrix, as `numpy.linalg.eig`. More...

class	einsum_error
	Thrown when an einsum specification is malformed or inconsistent. More...

class	hermitian_eigen
	The eigenvalues and eigenvectors of a complex Hermitian matrix. More...

struct	least_squares_solution
	A least-squares solution and what `numpy.linalg.lstsq` reports beside it (least_squares). More...

class	lu_decomposition
	An LU factorisation with partial pivoting of a square matrix. More...

class	matrix
	A dense, row-major matrix of `T` with value semantics. More...

class	non_finite_error
	Thrown when a decomposition's input holds an infinite or NaN entry: no factorisation of it means anything, and the iterative ones (eigh, SVD) would stop at once and return plausible values. More...

class	not_positive_definite_error
	Thrown when a Cholesky factorisation finds the matrix not positive-definite. More...

class	qr_decomposition
	A QR factorisation `A = Q R` of an `m × n` matrix, as NumPy's `qr`, reduced by default. More...

class	rank_cutoff
	The magnitude at or under which a singular value counts as zero, in a rank (matrix_rank), a pseudo-inverse (pinv) or a least-squares solution (lstsq): by default NumPy's, the largest singular value times the machine epsilon times the larger dimension; relative to the largest, as NumPy's `rcond` and `rtol`; or absolute, as `matrix_rank`'s `tol`. More...

class	singular_matrix_error
	Thrown when a matrix has no unique solution (it is singular). More...

class	svd_decomposition
	The singular value decomposition `A = U Σ Vᴴ` of an `m × n` matrix. More...

class	symmetric_eigen
	The eigenvalues and eigenvectors of a real symmetric matrix. More...

class	tensor
	A dense, row-major, N-dimensional tensor of `T` with value semantics. More...

class	vector
	A dense vector of `T` with value semantics. More...

Typedefs
template<typename T >
using	buffer = std::vector< T, detail::default_init_allocator< T > >
	The storage of vector, matrix and tensor, namely a std::vector whose `buffer<T>(n)` and `resize(n)` leave a `double` uninitialised, for a result about to be written in full. `buffer<T>(n, T {})` zeroes, as `std::vector<T>(n)` does.

Enumerations
enum class	qr_mode { reduced , complete }
	The factors a numerical::qr_decomposition keeps, as `numpy.linalg.qr`'s modes name them. More...

enum class	eigen_method { tridiagonal_ql , jacobi , divide_and_conquer }
	How numerical::symmetric_eigen diagonalises. More...

enum class	svd_method { bidiagonal_qr , jacobi , divide_and_conquer }
	How numerical::svd_decomposition computes the decomposition of `R`, after a Householder QR. More...

enum class	gemm_path : std::uint8_t { portable = 0 , neon = 1 , avx2 = 2 , avx512 = 3 , amx = 4 }
	The kernels a `double` matrix product can run on (a `float` one too on arm64), from the most portable. Each is used only where the processor has it: NEON on every arm64 chip, the AVX ones where `cpuid` reports them, the Apple AMX on an M1 to M4 chip whose self-test passes (see x86.hpp and amx.hpp). More...

Functions
template<typename T >
vector< T >	solve (matrix< T > a, const vector< T > &b)
	Solves the linear system `A x = b` for a square `A`.

template<typename T >
matrix< T >	solve (matrix< T > a, const matrix< T > &b)
	Solves `A X = B` for a square `A` and every column of `B` at once, as `numpy.linalg.solve` with a matrix right-hand side: one factorisation of `A`, then one substitution per column.

template<typename T >
T	determinant (matrix< T > a)
	The determinant of a square matrix.

template<typename T >
matrix< T >	inverse (matrix< T > a)
	The inverse of a square matrix.

template<typename T >
cholesky_decomposition< T >	cholesky (matrix< T > a)
	Cholesky-factorises a Hermitian positive-definite matrix.

template<typename T >
qr_decomposition< T >	qr (const matrix< T > &a)
	QR-factorises an `m × n` matrix.

template<typename T >
vector< T >	lstsq (const matrix< T > &a, const vector< T > &b, rank_cutoff cutoff={})
	The minimum-norm least-squares solution of `A x ≈ b`, as NumPy's `lstsq` gives it.

template<typename T >
matrix< T >	lstsq (const matrix< T > &a, const matrix< T > &b, rank_cutoff cutoff={})
	The minimum-norm least-squares solution of `A X ≈ B` for every column of `B` at once, as `numpy.linalg.lstsq` with a matrix right-hand side: one SVD of `A`, `X = V Σ⁺ Uᴴ B`; or, where a bound on the condition number shows the rank full with room to spare, the same solution to rounding by one QR factorisation (detail::lstsq_by_qr).

template<typename T >
vector< double >	singular_values (const matrix< T > &a)
	The singular values of `a`, descending, without its singular vectors, as `numpy.linalg.svd(a, compute_uv=False)`: the Householder `R` of a tall `A` (of `Aᴴ` when wide), its Golub–Kahan bidiagonal and the implicit-shift QR iteration on it, no basis formed nor rotated. Each value to an absolute accuracy of the largest times the machine epsilon, as svd_decomposition's.

template<typename T >
vector< double >	symmetric_eigenvalues (const matrix< T > &a)
	The eigenvalues of the real symmetric or complex Hermitian `a`, ascending, without its eigenvectors, as `numpy.linalg.eigvalsh`: its Householder tridiagonal (no basis formed) and the implicit QL iteration on it, no rotation kept. Only the lower triangle is read, and of the diagonal its real part. Each value to an absolute accuracy of the largest magnitude times the machine epsilon, as symmetric_eigen's.

template<typename T >
least_squares_solution< T >	least_squares (const matrix< T > &a, const matrix< T > &b, rank_cutoff cutoff={})
	lstsq, with the residuals, the rank and the singular values `numpy.linalg.lstsq` returns beside the solution: the residual sums of squares are given only where the solution is unique and need not be exact, `A` tall (`m > n`) and of full column rank, as NumPy gives them.

template<typename T >
std::size_t	matrix_rank (const matrix< T > &a, rank_cutoff cutoff={}, bool hermitian=false)
	The numerical rank of `a`, as `numpy.linalg.matrix_rank` counts it: its singular values above `cutoff`.

template<typename T >
matrix< T >	pinv (const matrix< T > &a, rank_cutoff cutoff={}, bool hermitian=false)
	The Moore–Penrose pseudo-inverse of the `m × n` `a` (`n × m`), as `numpy.linalg.pinv` gives it: `V Σ⁺ Uᴴ`, a singular value at or under `cutoff` counting as zero.

template<typename T >
vector< std::complex< double > >	eigenvalues (const matrix< T > &a)
	The eigenvalues alone of the square `a`, as `numpy.linalg.eigvals`.

gemm_path	limit_gemm_path (gemm_path most) noexcept
	Keeps `double` matrix products on `most` or a more portable path: the most capable path at or below it that can run here. `gemm_path::amx` restores the best one; `gemm_path::portable` keeps them on the compiler's own vectors, with no instruction chosen at run time. Paths below the cap are not probed.

gemm_path	gemm_path_in_use () noexcept
	The path `double` matrix products take.

template<typename T >
matrix< T >	operator+ (const matrix< T > &a, const matrix< T > &b)
	Element-wise sum.

template<typename T >
matrix< T >	operator- (const matrix< T > &a, const matrix< T > &b)
	Element-wise difference.

template<typename T >
matrix< T >	operator* (const T &scalar, const matrix< T > &m)
	Scalar times matrix.

template<typename T >
matrix< T >	operator* (const matrix< T > &m, const T &scalar)
	Matrix times scalar.

template<typename T >
vector< T >	operator* (const matrix< T > &a, const vector< T > &v)
	Matrix–vector product: each element its row's dot product with `v` (detail::matrix_times_vector).

template<typename T >
matrix< T >	operator* (const matrix< T > &a, const matrix< T > &b)
	Matrix–matrix product.

template<typename T >
matrix< T >	transpose (const matrix< T > &m)
	The transpose, by blocks of transpose_tile, the result's blocks of rows shared among the threads.

template<typename T >
bool	approx_equal (const matrix< T > &a, const matrix< T > &b, double tolerance=1e-9)
	Element-wise closeness within a tolerance (for inexact scalars).

void	set_max_threads (std::size_t threads) noexcept
	Caps the threads a kernel uses at `threads`, the caller's included (1: none beside it, and no thread is started); 0 lifts the cap. For a process that runs several kernels side by side (one per worker process, say), so that together they do not ask for more threads than the machine has.

std::size_t	max_threads ()
	The threads a kernel may use: the cap if one is set, else the hardware threads (starting none).

template<typename T >
bool	approx_equal (const tensor< T > &a, const tensor< T > &b, double tolerance=1e-9)
	Element-wise closeness within a tolerance (for inexact scalars).

template<typename T , typename ... Rest>
tensor< T >	einsum (const std::string &spec, const tensor< T > &first, const Rest &... rest)
	Einstein summation over one or more tensor operands.

template<typename T >
tensor< T >	einsum (const std::string &spec, const std::vector< const tensor< T > * > &operands)
	Einstein summation over operands counted at run time: einsum for a caller that holds them in a container (an interpreter evaluating `einsum("ij,jk->ik", a, b)` from a parse).

template<typename T >
tensor< T >	matmul (const tensor< T > &a, const tensor< T > &b)
	The matrix product of two tensors, as NumPy's `matmul`: each operand a stack of matrices along its leading axes, those axes broadcast against the other's (aligned from the right, an extent of 1 repeated), one product per matrix of the broadcast stack through the packed kernel, the products in parallel. A vector first is a row (its axis absent from the result), a vector second a column (likewise); two vectors give their dot product, of rank 0. The broadcast operands are not copied: a repeated matrix is read where it is.

std::vector< std::size_t >	broadcast_shape (const std::vector< std::size_t > &a, const std::vector< std::size_t > &b)
	The shape two shapes broadcast to, as NumPy's `broadcast_shapes`.

template<typename T , typename Op >
tensor< T >	broadcast_apply (const tensor< T > &a, const tensor< T > &b, const Op &op)
	`op(a, b)` element by element, the operands broadcast against each other.

template<detail::split_complex C>
tensor< C >	numpy_multiply (const tensor< C > &a, const tensor< C > &b)
	`a · b` of complex tensors, broadcast, each product as NumPy's array loop computes it (detail::numpy_multiply on this processor): the binding's ``, where `std::complex`'s operator (this library's ``) recovers infinities NumPy's makes NaN.

template<detail::split_complex C>
tensor< C >	numpy_divide (const tensor< C > &a, const tensor< C > &b)
	`a / b` of complex tensors, broadcast, each quotient as NumPy's array loop computes it (detail::numpy_divide on this processor): the binding's `/`.

template<typename T >
tensor< T >	operator+ (const tensor< T > &a, const tensor< T > &b)
	Element-wise sum, broadcast.

template<typename T >
tensor< T >	operator- (const tensor< T > &a, const tensor< T > &b)
	Element-wise difference, broadcast.

template<typename T >
tensor< T >	operator* (const tensor< T > &a, const tensor< T > &b)
	Element-wise product, broadcast.

template<typename T >
tensor< T >	operator/ (const tensor< T > &a, const tensor< T > &b)
	Element-wise quotient, broadcast.

template<typename T >
tensor< T >	operator- (const tensor< T > &a)
	Element-wise negation.

template<typename T , typename F >
auto	map (const tensor< T > &t, const F &f) -> tensor< std::decay_t< decltype(f(std::declval< const T & >()))> >
	`f(x)` for every element, in a tensor of the same shape and of `f`'s result type.

template<typename T >
tensor< T >	pow (const tensor< T > &a, const tensor< T > &b)
	`a` to the power `b` element by element, broadcast, as `std::pow`; a single exponent of 2, 0.5 or -1 is computed as NumPy computes it, `x·x`, `sqrt(x)` or `1/x` (exact where `pow` rounds: a complex square has no stray real part).

template<typename T >
T	sum (const tensor< T > &t)
	The sum of every element, by pairwise halving (0 for an empty tensor).

template<typename T >
tensor< T >	sum (const tensor< T > &t, std::size_t axis)
	The sums along `axis`, by pairwise halving.

template<typename T >
T	mean (const tensor< T > &t)
	The mean of every element. An empty tensor's mean is NaN, as NumPy's (which also warns).

template<typename T >
tensor< T >	mean (const tensor< T > &t, std::size_t axis)
	The means along `axis`.

template<std::totally_ordered T>
T	max (const tensor< T > &t)
	The largest element (a NaN, if there is one).

template<std::totally_ordered T>
T	min (const tensor< T > &t)
	The smallest element (a NaN, if there is one).

template<std::totally_ordered T>
tensor< T >	max (const tensor< T > &t, std::size_t axis)
	The largest elements along `axis`.

template<std::totally_ordered T>
tensor< T >	min (const tensor< T > &t, std::size_t axis)
	The smallest elements along `axis`.

template<typename T >
tensor< T >	reshape (const tensor< T > &t, std::vector< std::size_t > shape)
	The same elements in the same order under another shape.

template<typename T >
tensor< T >	full (std::vector< std::size_t > shape, const T &value)
	A tensor of `shape` with every element `value`.

template<typename T >
tensor< T >	zeros (std::vector< std::size_t > shape)
	A tensor of `shape` full of zeros.

template<typename T >
tensor< T >	ones (std::vector< std::size_t > shape)
	A tensor of `shape` full of ones.

template<typename T >
tensor< T >	eye (std::size_t n)
	The `n × n` identity.

tensor< double >	arange (double start, double stop, double step)
	`start, start + step, …` below `stop` (above it for a negative step), as NumPy's `arange`: its length is `ceil((stop − start) / step)` and element `i` is `start + i·delta` with `delta = (start + step) − start`, the step as the first two elements hold it, rounded once (a fused multiply-add, as NumPy's arm64 builds compute it; one without fusion rounds the product first and can differ in the last bits).

tensor< double >	linspace (double start, double stop, std::size_t count)
	`count` values evenly spaced from `start` to `stop`, both included, as NumPy's `linspace`: element `i` is `start + i·step` with `step = (stop − start) / (count − 1)`, and the last is `stop` exactly.

template<typename T >
vector< T >	operator+ (const vector< T > &a, const vector< T > &b)
	Element-wise sum.

template<typename T >
vector< T >	operator- (const vector< T > &a, const vector< T > &b)
	Element-wise difference.

template<typename T >
vector< T >	operator* (const T &scalar, const vector< T > &v)
	Scalar times vector.

template<typename T >
vector< T >	operator* (const vector< T > &v, const T &scalar)
	Vector times scalar.

template<typename T >
T	dot (const vector< T > &a, const vector< T > &b)
	Dot product.

template<typename T >
bool	approx_equal (const vector< T > &a, const vector< T > &b, double tolerance=1e-9)
	Element-wise closeness within a tolerance (for inexact scalars).

Variables
constexpr std::size_t	transpose_tile {32}
	The square blocks transpose copies one at a time: a block's rows of the source and of the result both stay in the L1 cache, where a row of the result walked a column of the source a stride apart.

Detailed Description

The NumericAL public API.

Typedef Documentation

◆ buffer

template<typename T >

using numerical::buffer = typedef std::vector<T, detail::default_init_allocator<T> >

The storage of vector, matrix and tensor, namely a std::vector whose buffer<T>(n) and resize(n) leave a double uninitialised, for a result about to be written in full. buffer<T>(n, T {}) zeroes, as std::vector<T>(n) does.

Definition at line 157 of file buffer.hpp.

Enumeration Type Documentation

◆ eigen_method

enum class numerical::eigen_method

strong

How numerical::symmetric_eigen diagonalises.

Enumerator
tridiagonal_ql	Householder tridiagonalisation then implicit QL, as LAPACK's `dsyev`: the vectors by rotations.
jacobi	Cyclic Jacobi rotations: several times slower, kept as an independent second method.
divide_and_conquer	Householder tridiagonalisation then divide and conquer, as LAPACK's `dsyevd`: the default.

Definition at line 4382 of file decomposition.hpp.

◆ gemm_path

enum class numerical::gemm_path : std::uint8_t

strong

The kernels a double matrix product can run on (a float one too on arm64), from the most portable. Each is used only where the processor has it: NEON on every arm64 chip, the AVX ones where cpuid reports them, the Apple AMX on an M1 to M4 chip whose self-test passes (see x86.hpp and amx.hpp).

Enumerator
portable	The compiler's vectors (16 bytes on x86-64 wheels), as the compiler rounds `c + a·b` (see gemm).
neon	arm64's vectors with fused multiply-adds, for `double` and `float`.
avx2	AVX2 with fused multiply-adds.
avx512	AVX-512 with fused multiply-adds.
amx	The Apple AMX coprocessor, fused multiply-adds; the products it does not suit on the NEON path's tile.

Definition at line 79 of file gemm.hpp.

◆ qr_mode

enum class numerical::qr_mode

strong

The factors a numerical::qr_decomposition keeps, as numpy.linalg.qr's modes name them.

Enumerator
reduced	`Q` is `m × k`, `R` `k × n`, `k = min(m, n)`: NumPy's default.
complete	`Q` is `m × m`, its columns past `k` a basis of the complement of `A`'s column space; `R` is `m × n`.

Definition at line 2535 of file decomposition.hpp.

◆ svd_method

enum class numerical::svd_method

strong

How numerical::svd_decomposition computes the decomposition of R, after a Householder QR.

Enumerator
bidiagonal_qr	Golub–Kahan bidiagonalisation then implicit-shift QR, as LAPACK's `dgesvd`.
jacobi	One-sided Jacobi: several times slower, each singular value to high relative accuracy.
divide_and_conquer	Golub–Kahan bidiagonalisation then divide and conquer, as LAPACK's `dgesdd`: the default.

Definition at line 5994 of file decomposition.hpp.

Function Documentation

◆ approx_equal() [1/3]

template<typename T >

bool numerical::approx_equal	(	const matrix< T > &	a,
		const matrix< T > &	b,
		double	tolerance = `1e-9`
	)

Element-wise closeness within a tolerance (for inexact scalars).

Parameters

[in]	a,b	Matrices to compare.
[in]	tolerance	Maximum allowed magnitude of each element-wise difference.

Returns: true if the shapes match and every element pair is equal or differs by at most tolerance; a NaN is close to nothing.

Definition at line 278 of file matrix.hpp.

◆ approx_equal() [2/3]

template<typename T >

bool numerical::approx_equal	(	const tensor< T > &	a,
		const tensor< T > &	b,
		double	tolerance = `1e-9`
	)

Element-wise closeness within a tolerance (for inexact scalars).

Parameters

[in]	a,b	Tensors to compare.
[in]	tolerance	Maximum allowed magnitude of each element-wise difference.

Returns: true if the shapes match and every element pair is equal or differs by at most tolerance; a NaN is close to nothing.

Definition at line 192 of file tensor.hpp.

◆ approx_equal() [3/3]

template<typename T >

bool numerical::approx_equal	(	const vector< T > &	a,
		const vector< T > &	b,
		double	tolerance = `1e-9`
	)

Element-wise closeness within a tolerance (for inexact scalars).

Parameters

[in]	a,b	Vectors to compare.
[in]	tolerance	Maximum allowed magnitude of each element-wise difference.

Returns: true if the sizes match and every element pair is equal or differs by at most tolerance; a NaN is close to nothing.

Definition at line 209 of file vector.hpp.

◆ arange()

tensor< double > numerical::arange	(	double	start,
		double	stop,
		double	step
	)

inline

start, start + step, … below stop (above it for a negative step), as NumPy's arange: its length is ceil((stop − start) / step) and element i is start + i·delta with delta = (start + step) − start, the step as the first two elements hold it, rounded once (a fused multiply-add, as NumPy's arm64 builds compute it; one without fusion rounds the product first and can differ in the last bits).

Exceptions

dimension_error If step is zero or not finite.

Definition at line 859 of file tensor_ops.hpp.

◆ broadcast_apply()

template<typename T , typename Op >

tensor< T > numerical::broadcast_apply	(	const tensor< T > &	a,
		const tensor< T > &	b,
		const Op &	op
	)

op(a, b) element by element, the operands broadcast against each other.

Exceptions

dimension_error If the shapes do not broadcast.

Definition at line 538 of file tensor_ops.hpp.

◆ broadcast_shape()

std::vector< std::size_t > numerical::broadcast_shape	(	const std::vector< std::size_t > &	a,
		const std::vector< std::size_t > &	b
	)

inline

The shape two shapes broadcast to, as NumPy's broadcast_shapes.

Exceptions

dimension_error If an axis has two extents, neither of them 1.

Definition at line 515 of file tensor_ops.hpp.

◆ cholesky()

template<typename T >

cholesky_decomposition< T > numerical::cholesky ( matrix< T > a )

Cholesky-factorises a Hermitian positive-definite matrix.

See also: cholesky_decomposition

Definition at line 1907 of file decomposition.hpp.

◆ determinant()

template<typename T >

T numerical::determinant ( matrix< T > a )

The determinant of a square matrix.

Exceptions

dimension_error If not square.

Definition at line 1533 of file decomposition.hpp.

◆ dot()

template<typename T >

T numerical::dot	(	const vector< T > &	a,
		const vector< T > &	b
	)

Dot product.

Exceptions

dimension_error If the sizes differ.

Definition at line 190 of file vector.hpp.

◆ eigenvalues()

template<typename T >

vector< std::complex< double > > numerical::eigenvalues ( const matrix< T > & a )

The eigenvalues alone of the square a, as numpy.linalg.eigvals.

Exceptions

dimension_error	If `a` is not square.
non_finite_error	If an entry is infinite or NaN.
convergence_error	If the QR iterations did not converge.

Definition at line 3159 of file eigen.hpp.

◆ einsum() [1/2]

template<typename T >

tensor< T > numerical::einsum	(	const std::string &	spec,
		const std::vector< const tensor< T > * > &	operands
	)

Einstein summation over operands counted at run time: einsum for a caller that holds them in a container (an interpreter evaluating einsum("ij,jk->ik", a, b) from a parse).

Parameters

[in]	spec	The einsum specification.
[in]	operands	The operands, in the spec's order (none null).

Returns: The resulting tensor.

Exceptions

einsum_error	As einsum, and if `operands` is empty.
dimension_error	As einsum.

Definition at line 1315 of file tensor.hpp.

◆ einsum() [2/2]

template<typename T , typename ... Rest>

tensor< T > numerical::einsum	(	const std::string &	spec,
		const tensor< T > &	first,
		const Rest &...	rest
	)

Einstein summation over one or more tensor operands.

Evaluates spec — input subscripts separated by ,, then ->, then the output subscripts — against the operands. Examples: "ij,jk->ik" (matrix product), "i,i->" (dot product), "ij->ji" (transpose), "ii->" (trace), "i,j->ij" (outer product), "ij->i" (row sums).

Parameters

[in]	spec	The einsum specification.
[in]	first	The first operand (fixes the scalar type).
[in]	rest	Any further operands.

Returns: The resulting tensor (rank-0 when the output subscripts are empty).

Exceptions

einsum_error	If the spec is malformed, the operand count or a rank disagrees with it, or an output label is invalid.
dimension_error	If a shared label labels axes of differing extents.

Definition at line 1297 of file tensor.hpp.

◆ eye()

template<typename T >

tensor< T > numerical::eye ( std::size_t n )

The n × n identity.

Definition at line 842 of file tensor_ops.hpp.

◆ full()

template<typename T >

tensor< T > numerical::full	(	std::vector< std::size_t >	shape,
		const T &	value
	)

A tensor of shape with every element value.

Definition at line 816 of file tensor_ops.hpp.

◆ gemm_path_in_use()

gemm_path numerical::gemm_path_in_use ( )

inlinenoexcept

The path double matrix products take.

Definition at line 3243 of file gemm.hpp.

◆ inverse()

template<typename T >

matrix< T > numerical::inverse ( matrix< T > a )

The inverse of a square matrix.

Parameters

[in] a The matrix.

Returns: Its inverse.

Exceptions

dimension_error	If `a` is not square.
singular_matrix_error	If `a` is singular.

Definition at line 1546 of file decomposition.hpp.

◆ least_squares()

template<typename T >

least_squares_solution< T > numerical::least_squares	(	const matrix< T > &	a,
		const matrix< T > &	b,
		rank_cutoff	cutoff = `{}`
	)

lstsq, with the residuals, the rank and the singular values numpy.linalg.lstsq returns beside the solution: the residual sums of squares are given only where the solution is unique and need not be exact, A tall (m > n) and of full column rank, as NumPy gives them.

Parameters

[in]	a	The `m × n` matrix.
[in]	b	The right-hand sides (`m × k`).
[in]	cutoff	Where a singular value counts as zero, in the solution and the rank.

Exceptions

dimension_error	If `b's` rows differ from `a's`.
non_finite_error	If an entry of `a` is infinite or NaN.

Definition at line 6607 of file decomposition.hpp.

◆ limit_gemm_path()

gemm_path numerical::limit_gemm_path ( gemm_path most )

inlinenoexcept

Keeps double matrix products on most or a more portable path: the most capable path at or below it that can run here. gemm_path::amx restores the best one; gemm_path::portable keeps them on the compiler's own vectors, with no instruction chosen at run time. Paths below the cap are not probed.

Returns: The path now in use.

Definition at line 3236 of file gemm.hpp.

◆ linspace()

tensor< double > numerical::linspace	(	double	start,
		double	stop,
		std::size_t	count
	)

inline

count values evenly spaced from start to stop, both included, as NumPy's linspace: element i is start + i·step with step = (stop − start) / (count − 1), and the last is stop exactly.

Definition at line 884 of file tensor_ops.hpp.

◆ lstsq() [1/2]

template<typename T >

matrix< T > numerical::lstsq	(	const matrix< T > &	a,
		const matrix< T > &	b,
		rank_cutoff	cutoff = `{}`
	)

The minimum-norm least-squares solution of A X ≈ B for every column of B at once, as numpy.linalg.lstsq with a matrix right-hand side: one SVD of A, X = V Σ⁺ Uᴴ B; or, where a bound on the condition number shows the rank full with room to spare, the same solution to rounding by one QR factorisation (detail::lstsq_by_qr).

Template Parameters

T	The scalar element type (`double` or `std::complex<double>`).

Parameters

[in]	a	The `m × n` matrix.
[in]	b	The right-hand sides (`m × k`).
[in]	cutoff	Where a singular value counts as zero.

Returns: X (n × k).

Exceptions

dimension_error	If `b`'s row count differs from `a`'s.
non_finite_error	If an entry of `a` is infinite or NaN.

Definition at line 6487 of file decomposition.hpp.

◆ lstsq() [2/2]

template<typename T >

vector< T > numerical::lstsq	(	const matrix< T > &	a,
		const vector< T > &	b,
		rank_cutoff	cutoff = `{}`
	)

The minimum-norm least-squares solution of A x ≈ b, as NumPy's lstsq gives it.

Template Parameters

T	The scalar element type (`double` or `std::complex<double>`).

From the SVD of A: x = V Σ⁺ Uᴴ b, where a singular value at or under cutoff (by default NumPy's, rank_cutoff) counts as zero. Any shape, any rank: an overdetermined full-rank system gives its least-squares solution, an underdetermined or rank-deficient one the solution of smallest norm.

Parameters

[in]	a	The `m × n` matrix.
[in]	b	The right-hand side (size `m`).
[in]	cutoff	Where a singular value counts as zero.

Returns: x (size n).

Exceptions

dimension_error	If `b`'s size differs from the row count.
non_finite_error	If an entry of `a` is infinite or NaN.

Definition at line 6465 of file decomposition.hpp.

◆ map()

template<typename T , typename F >

auto numerical::map	(	const tensor< T > &	t,
		const F &	f
	)		-> tensor<std::decay_t<decltype(f(std::declval<const T&>()))>>

f(x) for every element, in a tensor of the same shape and of f's result type.

Definition at line 689 of file tensor_ops.hpp.

◆ matmul()

template<typename T >

tensor< T > numerical::matmul	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

The matrix product of two tensors, as NumPy's matmul: each operand a stack of matrices along its leading axes, those axes broadcast against the other's (aligned from the right, an extent of 1 repeated), one product per matrix of the broadcast stack through the packed kernel, the products in parallel. A vector first is a row (its axis absent from the result), a vector second a column (likewise); two vectors give their dot product, of rank 0. The broadcast operands are not copied: a repeated matrix is read where it is.

Parameters

[in] a,b The operands, of rank 1 or more.

Returns: a @ b, of the broadcast leading axes, then the rows of a and the columns of b.

Exceptions

dimension_error If an operand is of rank 0, the inner extents differ, or the leading axes do not broadcast.

Definition at line 1332 of file tensor.hpp.

◆ matrix_rank()

template<typename T >

std::size_t numerical::matrix_rank	(	const matrix< T > &	a,
		rank_cutoff	cutoff = `{}`,
		bool	hermitian = `false`
	)

The numerical rank of a, as numpy.linalg.matrix_rank counts it: its singular values above cutoff.

Parameters

[in]	a	The matrix.
[in]	cutoff	Where a singular value counts as zero; by default NumPy's (rank_cutoff).
[in]	hermitian	Whether `a` is real symmetric or complex Hermitian, as NumPy's `hermitian=True`: its singular values are then the magnitudes of its eigenvalues, cheaper to find (in `double`, real or complex; any other scalar takes its singular values), and only its lower triangle is read.

Exceptions

dimension_error	If `hermitian` and `a` is not square.
non_finite_error	If an entry read is infinite or NaN.

Definition at line 6789 of file decomposition.hpp.

◆ max() [1/2]

template<std::totally_ordered T>

T numerical::max ( const tensor< T > & t )

The largest element (a NaN, if there is one).

Exceptions

dimension_error If t is empty.

Definition at line 767 of file tensor_ops.hpp.

◆ max() [2/2]

template<std::totally_ordered T>

tensor< T > numerical::max	(	const tensor< T > &	t,
		std::size_t	axis
	)

The largest elements along axis.

Exceptions

dimension_error If axis is not an axis or is empty.

Definition at line 781 of file tensor_ops.hpp.

◆ max_threads()

std::size_t numerical::max_threads ( )

inlineprivate

The threads a kernel may use: the cap if one is set, else the hardware threads (starting none).

Definition at line 732 of file parallel.hpp.

◆ mean() [1/2]

template<typename T >

T numerical::mean ( const tensor< T > & t )

The mean of every element. An empty tensor's mean is NaN, as NumPy's (which also warns).

Definition at line 746 of file tensor_ops.hpp.

◆ mean() [2/2]

template<typename T >

tensor< T > numerical::mean	(	const tensor< T > &	t,
		std::size_t	axis
	)

The means along axis.

Exceptions

dimension_error If axis is not an axis.

Definition at line 753 of file tensor_ops.hpp.

◆ min() [1/2]

template<std::totally_ordered T>

T numerical::min ( const tensor< T > & t )

The smallest element (a NaN, if there is one).

Exceptions

dimension_error If t is empty.

Definition at line 774 of file tensor_ops.hpp.

◆ min() [2/2]

template<std::totally_ordered T>

tensor< T > numerical::min	(	const tensor< T > &	t,
		std::size_t	axis
	)

The smallest elements along axis.

Exceptions

dimension_error If axis is not an axis or is empty.

Definition at line 789 of file tensor_ops.hpp.

◆ numpy_divide()

template<detail::split_complex C>

tensor< C > numerical::numpy_divide	(	const tensor< C > &	a,
		const tensor< C > &	b
	)

a / b of complex tensors, broadcast, each quotient as NumPy's array loop computes it (detail::numpy_divide on this processor): the binding's /.

Definition at line 638 of file tensor_ops.hpp.

◆ numpy_multiply()

template<detail::split_complex C>

tensor< C > numerical::numpy_multiply	(	const tensor< C > &	a,
		const tensor< C > &	b
	)

a · b of complex tensors, broadcast, each product as NumPy's array loop computes it (detail::numpy_multiply on this processor): the binding's *, where std::complex's operator (this library's *) recovers infinities NumPy's makes NaN.

Definition at line 628 of file tensor_ops.hpp.

◆ ones()

template<typename T >

tensor< T > numerical::ones ( std::vector< std::size_t > shape )

A tensor of shape full of ones.

Definition at line 835 of file tensor_ops.hpp.

◆ operator*() [1/7]

template<typename T >

matrix< T > numerical::operator*	(	const matrix< T > &	a,
		const matrix< T > &	b
	)

Matrix–matrix product.

Exceptions

dimension_error If shapes mismatch.

Definition at line 225 of file matrix.hpp.

◆ operator*() [2/7]

template<typename T >

vector< T > numerical::operator*	(	const matrix< T > &	a,
		const vector< T > &	v
	)

Matrix–vector product: each element its row's dot product with v (detail::matrix_times_vector).

Exceptions

dimension_error If shapes mismatch.

Definition at line 207 of file matrix.hpp.

◆ operator*() [3/7]

template<typename T >

matrix< T > numerical::operator*	(	const matrix< T > &	m,
		const T &	scalar
	)

Matrix times scalar.

Parameters

[in] m,scalar Operands.

Returns: m·scalar.

Definition at line 198 of file matrix.hpp.

◆ operator*() [4/7]

template<typename T >

matrix< T > numerical::operator*	(	const T &	scalar,
		const matrix< T > &	m
	)

Scalar times matrix.

Parameters

[in] scalar,m Operands.

Returns: scalar·m.

Definition at line 187 of file matrix.hpp.

◆ operator*() [5/7]

template<typename T >

vector< T > numerical::operator*	(	const T &	scalar,
		const vector< T > &	v
	)

Scalar times vector.

Parameters

[in] scalar,v Operands.

Returns: scalar·v.

Definition at line 171 of file vector.hpp.

◆ operator*() [6/7]

template<typename T >

tensor< T > numerical::operator*	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

Element-wise product, broadcast.

Exceptions

dimension_error If the shapes do not broadcast.

Definition at line 663 of file tensor_ops.hpp.

◆ operator*() [7/7]

template<typename T >

vector< T > numerical::operator*	(	const vector< T > &	v,
		const T &	scalar
	)

Vector times scalar.

Parameters

[in] v,scalar Operands.

Returns: v·scalar.

Definition at line 182 of file vector.hpp.

◆ operator+() [1/3]

template<typename T >

matrix< T > numerical::operator+	(	const matrix< T > &	a,
		const matrix< T > &	b
	)

Element-wise sum.

Exceptions

dimension_error If the shapes differ.

Definition at line 161 of file matrix.hpp.

◆ operator+() [2/3]

template<typename T >

tensor< T > numerical::operator+	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

Element-wise sum, broadcast.

Exceptions

dimension_error If the shapes do not broadcast.

Definition at line 647 of file tensor_ops.hpp.

◆ operator+() [3/3]

template<typename T >

vector< T > numerical::operator+	(	const vector< T > &	a,
		const vector< T > &	b
	)

Element-wise sum.

Exceptions

dimension_error If the sizes differ.

Definition at line 145 of file vector.hpp.

◆ operator-() [1/4]

template<typename T >

matrix< T > numerical::operator-	(	const matrix< T > &	a,
		const matrix< T > &	b
	)

Element-wise difference.

Exceptions

dimension_error If the shapes differ.

Definition at line 174 of file matrix.hpp.

◆ operator-() [2/4]

template<typename T >

tensor< T > numerical::operator- ( const tensor< T > & a )

Element-wise negation.

Definition at line 679 of file tensor_ops.hpp.

◆ operator-() [3/4]

template<typename T >

tensor< T > numerical::operator-	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

Element-wise difference, broadcast.

Exceptions

dimension_error If the shapes do not broadcast.

Definition at line 655 of file tensor_ops.hpp.

◆ operator-() [4/4]

template<typename T >

vector< T > numerical::operator-	(	const vector< T > &	a,
		const vector< T > &	b
	)

Element-wise difference.

Exceptions

dimension_error If the sizes differ.

Definition at line 158 of file vector.hpp.

◆ operator/()

template<typename T >

tensor< T > numerical::operator/	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

Element-wise quotient, broadcast.

Exceptions

dimension_error If the shapes do not broadcast.

Definition at line 671 of file tensor_ops.hpp.

◆ pinv()

template<typename T >

matrix< T > numerical::pinv	(	const matrix< T > &	a,
		rank_cutoff	cutoff = `{}`,
		bool	hermitian = `false`
	)

The Moore–Penrose pseudo-inverse of the m × n a (n × m), as numpy.linalg.pinv gives it: V Σ⁺ Uᴴ, a singular value at or under cutoff counting as zero.

Parameters

[in]	a	The matrix.
[in]	cutoff	Where a singular value counts as zero; by default NumPy's (rank_cutoff).
[in]	hermitian	Whether `a` is real symmetric or complex Hermitian, as NumPy's `hermitian=True`: from its eigendecomposition `V Λ Vᴴ`, `V Λ⁺ Vᴴ`, cheaper (in `double`, real or complex; any other scalar takes its singular value decomposition), and only its lower triangle is read.

Exceptions

dimension_error	If `hermitian` and `a` is not square.
non_finite_error	If an entry read is infinite or NaN.

Definition at line 6815 of file decomposition.hpp.

◆ pow()

template<typename T >

tensor< T > numerical::pow	(	const tensor< T > &	a,
		const tensor< T > &	b
	)

a to the power b element by element, broadcast, as std::pow; a single exponent of 2, 0.5 or -1 is computed as NumPy computes it, x·x, sqrt(x) or 1/x (exact where pow rounds: a complex square has no stray real part).

Definition at line 704 of file tensor_ops.hpp.

◆ qr()

template<typename T >

qr_decomposition< T > numerical::qr ( const matrix< T > & a )

QR-factorises an m × n matrix.

See also: qr_decomposition

Definition at line 2667 of file decomposition.hpp.

◆ reshape()

template<typename T >

tensor< T > numerical::reshape	(	const tensor< T > &	t,
		std::vector< std::size_t >	shape
	)

The same elements in the same order under another shape.

Exceptions

dimension_error If the shapes hold different counts of elements.

Definition at line 800 of file tensor_ops.hpp.

◆ set_max_threads()

void numerical::set_max_threads ( std::size_t threads )

inlineprivatenoexcept

Caps the threads a kernel uses at threads, the caller's included (1: none beside it, and no thread is started); 0 lifts the cap. For a process that runs several kernels side by side (one per worker process, say), so that together they do not ask for more threads than the machine has.

Definition at line 726 of file parallel.hpp.

◆ singular_values()

template<typename T >

vector< double > numerical::singular_values ( const matrix< T > & a )

The singular values of a, descending, without its singular vectors, as numpy.linalg.svd(a, compute_uv=False): the Householder R of a tall A (of Aᴴ when wide), its Golub–Kahan bidiagonal and the implicit-shift QR iteration on it, no basis formed nor rotated. Each value to an absolute accuracy of the largest times the machine epsilon, as svd_decomposition's.

Exceptions

non_finite_error	If an entry of `a` is infinite or NaN.
convergence_error	If the iteration does not converge within its bound (not observed).

Definition at line 6511 of file decomposition.hpp.

◆ solve() [1/2]

template<typename T >

matrix< T > numerical::solve	(	matrix< T >	a,
		const matrix< T > &	b
	)

Solves A X = B for a square A and every column of B at once, as numpy.linalg.solve with a matrix right-hand side: one factorisation of A, then one substitution per column.

Parameters

[in]	a	The square coefficient matrix (`n × n`).
[in]	b	The right-hand sides (`n × k`).

Returns: X (n × k).

Exceptions

dimension_error	If `a` is not square or `b`'s row count differs from it.
singular_matrix_error	If `a` is singular.

Definition at line 1525 of file decomposition.hpp.

◆ solve() [2/2]

template<typename T >

vector< T > numerical::solve	(	matrix< T >	a,
		const vector< T > &	b
	)

Solves the linear system A x = b for a square A.

Parameters

[in]	a	The square coefficient matrix.
[in]	b	The right-hand side.

Returns: The solution x.

Exceptions

dimension_error	If `a` is not square or `b`'s size differs from it.
singular_matrix_error	If `a` is singular.

Definition at line 1509 of file decomposition.hpp.

◆ sum() [1/2]

template<typename T >

T numerical::sum ( const tensor< T > & t )

The sum of every element, by pairwise halving (0 for an empty tensor).

Definition at line 726 of file tensor_ops.hpp.

◆ sum() [2/2]

template<typename T >

tensor< T > numerical::sum	(	const tensor< T > &	t,
		std::size_t	axis
	)

The sums along axis, by pairwise halving.

Exceptions

dimension_error If axis is not an axis.

Definition at line 733 of file tensor_ops.hpp.

◆ symmetric_eigenvalues()

template<typename T >

vector< double > numerical::symmetric_eigenvalues ( const matrix< T > & a )

The eigenvalues of the real symmetric or complex Hermitian a, ascending, without its eigenvectors, as numpy.linalg.eigvalsh: its Householder tridiagonal (no basis formed) and the implicit QL iteration on it, no rotation kept. Only the lower triangle is read, and of the diagonal its real part. Each value to an absolute accuracy of the largest magnitude times the machine epsilon, as symmetric_eigen's.

Exceptions

dimension_error	If `a` is not square.
non_finite_error	If an entry read is infinite or NaN.
convergence_error	If the iteration does not converge within its bound (not observed).

Definition at line 6549 of file decomposition.hpp.

◆ transpose()

template<typename T >

matrix< T > numerical::transpose ( const matrix< T > & m )

The transpose, by blocks of transpose_tile, the result's blocks of rows shared among the threads.

Parameters

[in] m The matrix.

Returns: An cols×rows matrix.

Definition at line 250 of file matrix.hpp.

◆ zeros()

template<typename T >

tensor< T > numerical::zeros ( std::vector< std::size_t > shape )

A tensor of shape full of zeros.

Definition at line 828 of file tensor_ops.hpp.

Variable Documentation

◆ transpose_tile

constexpr std::size_t numerical::transpose_tile {32}

inlineconstexpr

The square blocks transpose copies one at a time: a block's rows of the source and of the result both stay in the L1 cache, where a row of the result walked a column of the source a stride apart.

Definition at line 243 of file matrix.hpp.

Namespaces

Classes

Typedefs

Enumerations

Functions

Variables

Detailed Description

Typedef Documentation

◆ buffer

Enumeration Type Documentation

◆ eigen_method

◆ gemm_path

◆ qr_mode

◆ svd_method

Function Documentation

◆ approx_equal() [1/3]

◆ approx_equal() [2/3]

◆ approx_equal() [3/3]

◆ arange()

◆ broadcast_apply()

◆ broadcast_shape()

◆ cholesky()

◆ determinant()

◆ dot()

◆ eigenvalues()

◆ einsum() [1/2]

◆ einsum() [2/2]

◆ eye()

◆ full()

◆ gemm_path_in_use()

◆ inverse()

◆ least_squares()

◆ limit_gemm_path()

◆ linspace()

◆ lstsq() [1/2]

◆ lstsq() [2/2]

◆ map()

◆ matmul()

◆ matrix_rank()

◆ max() [1/2]

◆ max() [2/2]

◆ max_threads()

◆ mean() [1/2]

◆ mean() [2/2]

◆ min() [1/2]

◆ min() [2/2]

◆ numpy_divide()

◆ numpy_multiply()

◆ ones()

◆ operator*() [1/7]

◆ operator*() [2/7]

◆ operator*() [3/7]

◆ operator*() [4/7]

◆ operator*() [5/7]

◆ operator*() [6/7]

◆ operator*() [7/7]

◆ operator+() [1/3]

◆ operator+() [2/3]

◆ operator+() [3/3]

◆ operator-() [1/4]

◆ operator-() [2/4]

◆ operator-() [3/4]

◆ operator-() [4/4]

◆ operator/()

◆ pinv()

◆ pow()

◆ qr()

◆ reshape()

◆ set_max_threads()

◆ singular_values()

◆ solve() [1/2]

◆ solve() [2/2]

◆ sum() [1/2]

◆ sum() [2/2]

◆ symmetric_eigenvalues()

◆ transpose()

◆ zeros()

Variable Documentation

◆ transpose_tile