Namespaces
namespace	detail

Functions
template<typename T >
int	gemm (char transa, char transb, std::ptrdiff_t m, std::ptrdiff_t n, std::ptrdiff_t k, T alpha, const T a, std::ptrdiff_t lda, const T b, std::ptrdiff_t ldb, T beta, T *c, std::ptrdiff_t ldc)
	`C ← α op(A) op(B) + β C`, column-major, as BLAS's `dgemm`: `op(A)` is `m × k`, `op(B)` `k × n`, `C` `m × n`. With β zero, `C` is not read (a NaN there does not survive), as the reference's. Each term is `(α b) a`, α applied to B's element first, as the reference applies it; the sums run through NumericAL's product kernel, each element's terms in ascending order of the depth.

double	dot (std::ptrdiff_t n, const double x, std::ptrdiff_t incx, const double y, std::ptrdiff_t incy)
	`Σ x_i y_i` over `n` elements `incx` and `incy` apart (from the last where negative), as BLAS's `ddot`; 0 where `n ≤ 0`. Contiguous vectors by numerical::detail::chunked_dot (the vector units, the chunks in parallel, the same sum on any number of threads); others gathered first.

float	dot (std::ptrdiff_t n, const float x, std::ptrdiff_t incx, const float y, std::ptrdiff_t incy)
	dot of `float` elements, summed in `float`, as BLAS's `sdot`.

template<typename T >
std::complex< T >	complex_dot (bool conjugate, std::ptrdiff_t n, const std::complex< T > x, std::ptrdiff_t incx, const std::complex< T > y, std::ptrdiff_t incy) noexcept
	`Σ x_i y_i` (`conjugate` false, BLAS's `zdotu`) or `Σ conj(x_i) y_i` (true, `zdotc`) over `n` complex elements `incx` and `incy` apart (from the last where negative); 0 where `n ≤ 0`. Contiguous vectors (both apart by 1) go by vectors of their parts, as the library's own complex dot products (numerical::detail::complex_dot_product, numerical::detail::conjugate_dot_product); others term by term in the order of their index. Either on the calling thread.

int	gemv (char trans, std::ptrdiff_t m, std::ptrdiff_t n, double alpha, const double a, std::ptrdiff_t lda, const double x, std::ptrdiff_t incx, double beta, double *y, std::ptrdiff_t incy)
	`y ← α op(A) x + β y`, column-major `A` `m × n`, as BLAS's `dgemv`: the vectors `incx` and `incy` apart (from the last element where negative), `y` not read where `β = 0`. `Aᵀ x` is each column's dot product with `x` (numerical::detail::row_dots, `A`'s columns the rows of `Aᵀ`); `A x` adds `α x_j` times each column, by groups of gemv_rows rows in parallel.

int	syrk (char uplo, char trans, std::ptrdiff_t n, std::ptrdiff_t k, double alpha, const double a, std::ptrdiff_t lda, double beta, double c, std::ptrdiff_t ldc)
	`C ← α op(A) op(A)ᵀ + β C` on the triangle `uplo` (`U` or `L`) of the column-major `n × n` `C`, `op(A)` `n × k`, as BLAS's `dsyrk` (NumPy's `a @ a.T`): the other triangle neither read nor written, `C` not read where `β = 0`. By blocks of syrk_block columns: a diagonal block formed whole aside and its triangle kept, the rows past it (`L`) or before it (`U`) by gemm.

int	trsm (char side, char uplo, char transa, char diag, std::ptrdiff_t m, std::ptrdiff_t n, double alpha, const double a, std::ptrdiff_t lda, double b, std::ptrdiff_t ldb)
	`B ← α op(A)⁻¹ B` (`side` `L`) or `α B op(A)⁻¹` (`R`), column-major, `B` `m × n` and the triangular `A` (`uplo` `U` or `L`, its other triangle not read; a `diag` `U`nit diagonal not read either) `m × m` or `n × n`, as BLAS's `dtrsm`; `A` not read where `α = 0`. Every case is ‘A’ X = B'`with`A'`triangular:` A`or a transposed copy of it,`B`or a transposed copy of it (`X op(A) = B`is`op(A)ᵀ Xᵀ = Bᵀ`), solved by \ref detail::triangular_solve. A zero on the diagonal gives infinities, as the reference's division does. \return 0, or the position of the first illegal argument, as the reference passes it to`XERBLA`.

template<typename T , typename Index >
int	getrs (char trans, std::ptrdiff_t n, std::ptrdiff_t nrhs, const T a, std::ptrdiff_t lda, const Index ipiv, T *b, std::ptrdiff_t ldb)
	`A X = B`, `Aᵀ X = B` or `Aᴴ X = B` (`trans` `N`, `T` or `C`), column-major, from the factors and interchanges of a `dgetrf` (by any implementation), as LAPACK's `dgetrs`: `B` (`n × nrhs`) overwritten by `X`, each column by forward then backward substitution.

template<typename T , typename Index >
int	gesv (std::ptrdiff_t n, std::ptrdiff_t nrhs, T a, std::ptrdiff_t lda, Index ipiv, T *b, std::ptrdiff_t ldb)
	`A X = B`, column-major, as LAPACK's `dgesv`: `A` (`n × n`, `lda` apart) is factorised as `P A = L U` by NumericAL's blocked LU with partial pivoting, the factors left in `a` as `dgetrf` leaves them (`U` on and above the diagonal, `L` below, its unit diagonal not stored) and the row interchanges in `ipiv` (1-based: row `i` was interchanged with row `ipiv[i]`, in order); then `B` (`n × nrhs`, `ldb` apart) is overwritten by `X`. A zero pivot leaves `B` as it was. A NaN or an infinity in `A` or `B`, which NumericAL's LU refuses, is carried through as the reference's routines carry it: `A` factorised column by column (`dgetf2`), then `B` solved by getrs.

template<typename T >
int	potrf (char uplo, std::ptrdiff_t n, T *a, std::ptrdiff_t lda)
	The Cholesky factor of a symmetric positive-definite `A`, column-major, as LAPACK's `dpotrf`: with `uplo` `L`, `A = L Lᵀ` read from and written to `A`'s lower triangle; with `U`, `A = Uᵀ U` and the upper one. The other triangle is not referenced. Where `A` is not positive definite, `A` is left as it was (LAPACK leaves the factorisation it reached). A NaN or an infinity in the triangle read, which NumericAL's Cholesky refuses, is carried through as the reference's `dpotf2` carries it (detail::column_cholesky), its `INFO` too.

template<typename T , typename Index >
int	getrf (std::ptrdiff_t m, std::ptrdiff_t n, T a, std::ptrdiff_t lda, Index ipiv)
	`P A = L U` of the column-major `m × n` `A` in place, as LAPACK's `dgetrf`: the factors as `dgetrf` leaves them and the row interchanges in `ipiv` (1-based, `min(m, n)` of them). A square `A` by NumericAL's blocked LU, a rectangular one, or one holding a NaN or an infinity, column by column as the reference's `dgetf2` (which carries the NaN through).

template<typename T >
int	potrs (char uplo, std::ptrdiff_t n, std::ptrdiff_t nrhs, const T a, std::ptrdiff_t lda, T b, std::ptrdiff_t ldb)
	`A X = B` from the Cholesky factor of a `dpotrf` (by any implementation), column-major, as LAPACK's `dpotrs`: `uplo` `L` for `A = L Lᵀ` in the lower triangle, `U` for `A = Uᵀ U` in the upper one; `B` (`n × nrhs`) overwritten by `X`.

template<typename Index >
int	syevd (char jobz, char uplo, std::ptrdiff_t n, double a, std::ptrdiff_t lda, double w, double work, std::ptrdiff_t lwork, Index iwork, std::ptrdiff_t liwork)
	The eigenvalues, and with `jobz` `V` the eigenvectors, of a symmetric column-major `A`, as LAPACK's `dsyevd`: the triangle `uplo` (`L` or `U`) read; the eigenvalues ascending in `w`; the orthonormal eigenvectors overwriting `A`, column `j` that of `w[j]`; with `N`, `A` left as it was (LAPACK leaves it overwritten). NumericAL's divide and conquer computes them; the workspaces are not used, but their sizes are `dsyevd`'s: an `lwork` or `liwork` of -1 asks them (written to `work[0]` and `iwork[0]`), and a smaller one is illegal. A NaN or an infinity in the triangle read gives NaN eigenvalues (and eigenvectors) and `INFO` `n`, an eigenvalue that did not converge, as the reference's gives for a NaN (an infinity it turns to NaN eigenvalues, `INFO` 0 or positive).

template<typename Index >
int	gesdd (char jobz, std::ptrdiff_t m, std::ptrdiff_t n, double a, std::ptrdiff_t lda, double s, double u, std::ptrdiff_t ldu, double vt, std::ptrdiff_t ldvt, double work, std::ptrdiff_t lwork, Index iwork)
	The singular value decomposition `A = U Σ Vᵀ` of a column-major `m × n` `A`, as LAPACK's `dgesdd`: the singular values descending in `s`; with `jobz` `A`, all of `U` (`m × m`) and `Vᵀ` (`n × n`); with `S`, their first `min(m, n)` columns and rows; with `O`, those overwriting `A` (`U`'s where `m ≥ n`, else `Vᵀ`'s) and the other in `u` or `vt`; with `N`, the values alone (`A` then left as it was). NumericAL's divide and conquer computes them; the workspaces are not used, but their sizes are `dgesdd`'s: an `lwork` of -1 asks it (written to `work[0]`), and a smaller one is illegal.

int	geqrf (std::ptrdiff_t m, std::ptrdiff_t n, double a, std::ptrdiff_t lda, double tau, double *work, std::ptrdiff_t lwork)
	The QR factorisation `A = Q R` of a column-major `m × n` `A`, as LAPACK's `dgeqrf`: `R` overwrites the upper trapezoid of `A`; `Q = H₀ H₁ … H_{k-1}` (`k = min(m, n)`) is left as its reflectors `H_i = I − τ_i w wᵀ`, `w` 1 at row `i` and below it held in `A`'s column `i`, `τ_i` in `tau` (0 for the identity, which a column with nothing below the diagonal takes, as in `dlarfg`) — the form `dorgqr` and `dormqr` read. NumericAL's blocked Householder reflections compute them; the workspace is not used, but its size is `dgeqrf`'s least, `n`: an `lwork` of -1 asks it (written to `work[0]`), and a smaller one is illegal. A NaN or an infinity in `A` is carried through, as the reference's `dgeqrf` carries it.

int	orgqr (std::ptrdiff_t m, std::ptrdiff_t n, std::ptrdiff_t k, double a, std::ptrdiff_t lda, const double tau, double *work, std::ptrdiff_t lwork)
	The first `n` columns of `Q = H₀ H₁ … H_{k-1}`, from the reflectors of a `dgeqrf` (geqrf's or any implementation's: column `i` of `A` below row `i`, `τ_i` in `tau`), overwriting the column-major `m × n` `A`, as LAPACK's `dorgqr`: the reflectors gathered by blocks into `I − V T Vᵀ` (`T` as `dlarft` forms it) and applied through the matrix product kernel. The workspace is not used, but its size is `dorgqr`'s least, `n`: an `lwork` of -1 asks it (written to `work[0]`), and a smaller one is illegal. A NaN or an infinity in a reflector or in `tau` is carried through, as the reference's `dorgqr` carries it.

template<typename Index >
int	gelsd (std::ptrdiff_t m, std::ptrdiff_t n, std::ptrdiff_t nrhs, const double a, std::ptrdiff_t lda, double b, std::ptrdiff_t ldb, double s, double rcond, Index rank, double work, std::ptrdiff_t lwork, Index iwork)
	The minimum-norm least-squares solutions of `A X ≈ B`, column-major `A` `m × n` and `B` `max(m, n) × nrhs`, as LAPACK's `dgelsd`: a singular value at or under `rcond` times the largest counts as zero (`rcond ≤ 0` or `≥ 1`: lapack_precision); `X` overwrites `B`'s first `n` rows, the singular values (descending) go to `s` and the effective rank to `rank`. Where `m > n`, `B`'s rows `n…m-1` of each column hold its residual as LAPACK leaves it, a vector whose norm is `‖b − A x‖` (here that norm, then zeros), from which NumPy's `lstsq` takes its residuals. `A` is left as it was. NumericAL's least_squares computes them; the workspaces are not used, but their sizes are `dgelsd`'s: an `lwork` of -1 asks them (written to `work[0]` and `iwork[0]`), and a smaller `lwork` is illegal.

int	geev (char jobvl, char jobvr, std::ptrdiff_t n, const double a, std::ptrdiff_t lda, double wr, double wi, double vl, std::ptrdiff_t ldvl, double vr, std::ptrdiff_t ldvr, double work, std::ptrdiff_t lwork)
	The eigenvalues, and with `jobvl` or `jobvr` `V` the left or right eigenvectors, of a column-major `n × n` `A`, as LAPACK's `dgeev`: the eigenvalues as `wr` + i `wi`, a conjugate pair adjacent with its positive imaginary part first; a right eigenvector `v` (`A v = λ v`) and a left one `u` (`uᴴ A = λ uᴴ`) of unit norm with its largest component real, in `VR` and `VL` as `dgeev` packs them (a real eigenvalue's column its vector, a pair's two columns the real and imaginary parts of the first's). NumericAL's eigen_decomposition computes them; the left vectors are `Aᵀ`'s right ones, conjugated, each matched to the nearest unmatched eigenvalue. `A` is left as it was. The workspace is not used, but its size is `dgeev`'s least, `3n` (`4n` with vectors): an `lwork` of -1 asks it (written to `work[0]`), and a smaller one is illegal.

Variables
constexpr double	lapack_precision {0x1p-53}
	LAPACK's relative machine precision, ‘dlamch('E’)`: the unit roundoff`2⁻⁵³`, half of` std::numeric_limits<double>::epsilon()`. The default rank cutoff of`dgelsd` is relative to it.

constexpr std::size_t	syrk_block {128}
	The columns of syrk's smallest diagonal blocks, formed whole: their wasted half, at most `64 n k` multiply-adds, stays small beside the `n² k / 2` of the triangle.

constexpr std::size_t	syrk_kept_limit {(std::size_t {32} << 20U) / sizeof(double)}
	The most elements of syrk's operand copy a thread keeps between calls (32 MiB).

constexpr std::size_t	gemv_rows {64}
	The rows (or, for `Aᵀ x`, the columns) each task of gemv takes, 64 so that an order of a few hundred shares out (256 left order 300 two tasks).

constexpr std::size_t	gemv_weight {4}
	The weight of gemv's multiply-adds against a compute-bound product's: each reads an element of `A` once, so the product is bound by the memory, and spreads over more threads' caches for its work (six x86-64 AVX2 cores, `Aᵀ x`: order 300 17.7 → 7.5 µs, order 1000 74 → 60 µs).

constexpr std::size_t	trsm_block {64}
	The rows of trsm's blocks: within one, each column solved alone; past it, the matrix product kernel.

Function Documentation

◆ complex_dot()

template<typename T >

std::complex< T > numerical::blas::complex_dot	(	bool	conjugate,
		std::ptrdiff_t	n,
		const std::complex< T > *	x,
		std::ptrdiff_t	incx,
		const std::complex< T > *	y,
		std::ptrdiff_t	incy
	)

noexcept

Σ x_i y_i (conjugate false, BLAS's zdotu) or Σ conj(x_i) y_i (true, zdotc) over n complex elements incx and incy apart (from the last where negative); 0 where n ≤ 0. Contiguous vectors (both apart by 1) go by vectors of their parts, as the library's own complex dot products (numerical::detail::complex_dot_product, numerical::detail::conjugate_dot_product); others term by term in the order of their index. Either on the calling thread.

Definition at line 600 of file blas.hpp.

◆ dot() [1/2]

double numerical::blas::dot	(	std::ptrdiff_t	n,
		const double *	x,
		std::ptrdiff_t	incx,
		const double *	y,
		std::ptrdiff_t	incy
	)

inline

Σ x_i y_i over n elements incx and incy apart (from the last where negative), as BLAS's ddot; 0 where n ≤ 0. Contiguous vectors by numerical::detail::chunked_dot (the vector units, the chunks in parallel, the same sum on any number of threads); others gathered first.

Definition at line 573 of file blas.hpp.

◆ dot() [2/2]

float numerical::blas::dot	(	std::ptrdiff_t	n,
		const float *	x,
		std::ptrdiff_t	incx,
		const float *	y,
		std::ptrdiff_t	incy
	)

inline

dot of float elements, summed in float, as BLAS's sdot.

Definition at line 583 of file blas.hpp.

◆ geev()

int numerical::blas::geev	(	char	jobvl,
		char	jobvr,
		std::ptrdiff_t	n,
		const double *	a,
		std::ptrdiff_t	lda,
		double *	wr,
		double *	wi,
		double *	vl,
		std::ptrdiff_t	ldvl,
		double *	vr,
		std::ptrdiff_t	ldvr,
		double *	work,
		std::ptrdiff_t	lwork
	)

inline

The eigenvalues, and with jobvl or jobvr V the left or right eigenvectors, of a column-major n × n A, as LAPACK's dgeev: the eigenvalues as wr + i wi, a conjugate pair adjacent with its positive imaginary part first; a right eigenvector v (A v = λ v) and a left one u (uᴴ A = λ uᴴ) of unit norm with its largest component real, in VR and VL as dgeev packs them (a real eigenvalue's column its vector, a pair's two columns the real and imaginary parts of the first's). NumericAL's eigen_decomposition computes them; the left vectors are Aᵀ's right ones, conjugated, each matched to the nearest unmatched eigenvalue. A is left as it was. The workspace is not used, but its size is dgeev's least, 3n (4n with vectors): an lwork of -1 asks it (written to work[0]), and a smaller one is illegal.

Returns: INFO: 0, or -i where argument i is illegal (-4: a NaN or an infinity in A).

Definition at line 1922 of file blas.hpp.

◆ gelsd()

template<typename Index >

int numerical::blas::gelsd	(	std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		std::ptrdiff_t	nrhs,
		const double *	a,
		std::ptrdiff_t	lda,
		double *	b,
		std::ptrdiff_t	ldb,
		double *	s,
		double	rcond,
		Index *	rank,
		double *	work,
		std::ptrdiff_t	lwork,
		Index *	iwork
	)

The minimum-norm least-squares solutions of A X ≈ B, column-major A m × n and B max(m, n) × nrhs, as LAPACK's dgelsd: a singular value at or under rcond times the largest counts as zero (rcond ≤ 0 or ≥ 1: lapack_precision); X overwrites B's first n rows, the singular values (descending) go to s and the effective rank to rank. Where m > n, B's rows n…m-1 of each column hold its residual as LAPACK leaves it, a vector whose norm is ‖b − A x‖ (here that norm, then zeros), from which NumPy's lstsq takes its residuals. A is left as it was. NumericAL's least_squares computes them; the workspaces are not used, but their sizes are dgelsd's: an lwork of -1 asks them (written to work[0] and iwork[0]), and a smaller lwork is illegal.

Returns: INFO: 0, or -i where argument i is illegal (-4, -6: a NaN or an infinity in A, in B).

Definition at line 1797 of file blas.hpp.

◆ gemm()

template<typename T >

int numerical::blas::gemm	(	char	transa,
		char	transb,
		std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		std::ptrdiff_t	k,
		T	alpha,
		const T *	a,
		std::ptrdiff_t	lda,
		const T *	b,
		std::ptrdiff_t	ldb,
		T	beta,
		T *	c,
		std::ptrdiff_t	ldc
	)

C ← α op(A) op(B) + β C, column-major, as BLAS's dgemm: op(A) is m × k, op(B) k × n, C m × n. With β zero, C is not read (a NaN there does not survive), as the reference's. Each term is (α b) a, α applied to B's element first, as the reference applies it; the sums run through NumericAL's product kernel, each element's terms in ascending order of the depth.

Returns: 0, or the position of the first illegal argument as XERBLA numbers it (then nothing is computed).

Definition at line 468 of file blas.hpp.

◆ gemv()

int numerical::blas::gemv	(	char	trans,
		std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		double	alpha,
		const double *	a,
		std::ptrdiff_t	lda,
		const double *	x,
		std::ptrdiff_t	incx,
		double	beta,
		double *	y,
		std::ptrdiff_t	incy
	)

inline

y ← α op(A) x + β y, column-major A m × n, as BLAS's dgemv: the vectors incx and incy apart (from the last element where negative), y not read where β = 0. Aᵀ x is each column's dot product with x (numerical::detail::row_dots, A's columns the rows of Aᵀ); A x adds α x_j times each column, by groups of gemv_rows rows in parallel.

Returns: 0, or the position of the first illegal argument, as the reference passes it to XERBLA.

Definition at line 638 of file blas.hpp.

◆ geqrf()

int numerical::blas::geqrf	(	std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		double *	a,
		std::ptrdiff_t	lda,
		double *	tau,
		double *	work,
		std::ptrdiff_t	lwork
	)

inline

The QR factorisation A = Q R of a column-major m × n A, as LAPACK's dgeqrf: R overwrites the upper trapezoid of A; Q = H₀ H₁ … H_{k-1} (k = min(m, n)) is left as its reflectors H_i = I − τ_i w wᵀ, w 1 at row i and below it held in A's column i, τ_i in tau (0 for the identity, which a column with nothing below the diagonal takes, as in dlarfg) — the form dorgqr and dormqr read. NumericAL's blocked Householder reflections compute them; the workspace is not used, but its size is dgeqrf's least, n: an lwork of -1 asks it (written to work[0]), and a smaller one is illegal. A NaN or an infinity in A is carried through, as the reference's dgeqrf carries it.

Returns: INFO: 0, or -i where argument i is illegal.

Definition at line 1658 of file blas.hpp.

◆ gesdd()

template<typename Index >

int numerical::blas::gesdd	(	char	jobz,
		std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		double *	a,
		std::ptrdiff_t	lda,
		double *	s,
		double *	u,
		std::ptrdiff_t	ldu,
		double *	vt,
		std::ptrdiff_t	ldvt,
		double *	work,
		std::ptrdiff_t	lwork,
		Index *	iwork
	)

The singular value decomposition A = U Σ Vᵀ of a column-major m × n A, as LAPACK's dgesdd: the singular values descending in s; with jobz A, all of U (m × m) and Vᵀ (n × n); with S, their first min(m, n) columns and rows; with O, those overwriting A (U's where m ≥ n, else Vᵀ's) and the other in u or vt; with N, the values alone (A then left as it was). NumericAL's divide and conquer computes them; the workspaces are not used, but their sizes are dgesdd's: an lwork of -1 asks it (written to work[0]), and a smaller one is illegal.

Returns: INFO: 0, or -i where argument i is illegal (-4: a NaN or an infinity in A).

Definition at line 1547 of file blas.hpp.

◆ gesv()

template<typename T , typename Index >

int numerical::blas::gesv	(	std::ptrdiff_t	n,
		std::ptrdiff_t	nrhs,
		T *	a,
		std::ptrdiff_t	lda,
		Index *	ipiv,
		T *	b,
		std::ptrdiff_t	ldb
	)

A X = B, column-major, as LAPACK's dgesv: A (n × n, lda apart) is factorised as P A = L U by NumericAL's blocked LU with partial pivoting, the factors left in a as dgetrf leaves them (U on and above the diagonal, L below, its unit diagonal not stored) and the row interchanges in ipiv (1-based: row i was interchanged with row ipiv[i], in order); then B (n × nrhs, ldb apart) is overwritten by X. A zero pivot leaves B as it was. A NaN or an infinity in A or B, which NumericAL's LU refuses, is carried through as the reference's routines carry it: A factorised column by column (dgetf2), then B solved by getrs.

Returns: INFO: 0; i > 0 where U(i, i) is exactly zero; -i where argument i is illegal.

Definition at line 1077 of file blas.hpp.

◆ getrf()

template<typename T , typename Index >

int numerical::blas::getrf	(	std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		T *	a,
		std::ptrdiff_t	lda,
		Index *	ipiv
	)

P A = L U of the column-major m × n A in place, as LAPACK's dgetrf: the factors as dgetrf leaves them and the row interchanges in ipiv (1-based, min(m, n) of them). A square A by NumericAL's blocked LU, a rectangular one, or one holding a NaN or an infinity, column by column as the reference's dgetf2 (which carries the NaN through).

Returns: INFO: 0; i > 0 where U(i, i) is exactly zero; -i where argument i is illegal.

Definition at line 1267 of file blas.hpp.

◆ getrs()

template<typename T , typename Index >

int numerical::blas::getrs	(	char	trans,
		std::ptrdiff_t	n,
		std::ptrdiff_t	nrhs,
		const T *	a,
		std::ptrdiff_t	lda,
		const Index *	ipiv,
		T *	b,
		std::ptrdiff_t	ldb
	)

A X = B, Aᵀ X = B or Aᴴ X = B (trans N, T or C), column-major, from the factors and interchanges of a dgetrf (by any implementation), as LAPACK's dgetrs: B (n × nrhs) overwritten by X, each column by forward then backward substitution.

Returns: INFO: 0, or -i where argument i is illegal.

Definition at line 1320 of file blas.hpp.

◆ orgqr()

int numerical::blas::orgqr	(	std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		std::ptrdiff_t	k,
		double *	a,
		std::ptrdiff_t	lda,
		const double *	tau,
		double *	work,
		std::ptrdiff_t	lwork
	)

inline

The first n columns of Q = H₀ H₁ … H_{k-1}, from the reflectors of a dgeqrf (geqrf's or any implementation's: column i of A below row i, τ_i in tau), overwriting the column-major m × n A, as LAPACK's dorgqr: the reflectors gathered by blocks into I − V T Vᵀ (T as dlarft forms it) and applied through the matrix product kernel. The workspace is not used, but its size is dorgqr's least, n: an lwork of -1 asks it (written to work[0]), and a smaller one is illegal. A NaN or an infinity in a reflector or in tau is carried through, as the reference's dorgqr carries it.

Returns: INFO: 0, or -i where argument i is illegal.

Definition at line 1728 of file blas.hpp.

◆ potrf()

template<typename T >

int numerical::blas::potrf	(	char	uplo,
		std::ptrdiff_t	n,
		T *	a,
		std::ptrdiff_t	lda
	)

The Cholesky factor of a symmetric positive-definite A, column-major, as LAPACK's dpotrf: with uplo L, A = L Lᵀ read from and written to A's lower triangle; with U, A = Uᵀ U and the upper one. The other triangle is not referenced. Where A is not positive definite, A is left as it was (LAPACK leaves the factorisation it reached). A NaN or an infinity in the triangle read, which NumericAL's Cholesky refuses, is carried through as the reference's dpotf2 carries it (detail::column_cholesky), its INFO too.

Returns: INFO: 0; i > 0 where the leading minor of order i is not positive definite (or holds a NaN); -i where argument i is illegal.

Definition at line 1201 of file blas.hpp.

◆ potrs()

template<typename T >

int numerical::blas::potrs	(	char	uplo,
		std::ptrdiff_t	n,
		std::ptrdiff_t	nrhs,
		const T *	a,
		std::ptrdiff_t	lda,
		T *	b,
		std::ptrdiff_t	ldb
	)

A X = B from the Cholesky factor of a dpotrf (by any implementation), column-major, as LAPACK's dpotrs: uplo L for A = L Lᵀ in the lower triangle, U for A = Uᵀ U in the upper one; B (n × nrhs) overwritten by X.

Returns: INFO: 0, or -i where argument i is illegal.

Definition at line 1397 of file blas.hpp.

◆ syevd()

template<typename Index >

int numerical::blas::syevd	(	char	jobz,
		char	uplo,
		std::ptrdiff_t	n,
		double *	a,
		std::ptrdiff_t	lda,
		double *	w,
		double *	work,
		std::ptrdiff_t	lwork,
		Index *	iwork,
		std::ptrdiff_t	liwork
	)

The eigenvalues, and with jobz V the eigenvectors, of a symmetric column-major A, as LAPACK's dsyevd: the triangle uplo (L or U) read; the eigenvalues ascending in w; the orthonormal eigenvectors overwriting A, column j that of w[j]; with N, A left as it was (LAPACK leaves it overwritten). NumericAL's divide and conquer computes them; the workspaces are not used, but their sizes are dsyevd's: an lwork or liwork of -1 asks them (written to work[0] and iwork[0]), and a smaller one is illegal. A NaN or an infinity in the triangle read gives NaN eigenvalues (and eigenvectors) and INFO n, an eigenvalue that did not converge, as the reference's gives for a NaN (an infinity it turns to NaN eigenvalues, INFO 0 or positive).

Returns: INFO: 0; n where A holds a NaN or an infinity; -i where argument i is illegal.

Definition at line 1457 of file blas.hpp.

◆ syrk()

int numerical::blas::syrk	(	char	uplo,
		char	trans,
		std::ptrdiff_t	n,
		std::ptrdiff_t	k,
		double	alpha,
		const double *	a,
		std::ptrdiff_t	lda,
		double	beta,
		double *	c,
		std::ptrdiff_t	ldc
	)

inline

C ← α op(A) op(A)ᵀ + β C on the triangle uplo (U or L) of the column-major n × n C, op(A) n × k, as BLAS's dsyrk (NumPy's a @ a.T): the other triangle neither read nor written, C not read where β = 0. By blocks of syrk_block columns: a diagonal block formed whole aside and its triangle kept, the rows past it (L) or before it (U) by gemm.

Returns: 0, or the position of the first illegal argument, as the reference passes it to XERBLA.

Definition at line 828 of file blas.hpp.

◆ trsm()

int numerical::blas::trsm	(	char	side,
		char	uplo,
		char	transa,
		char	diag,
		std::ptrdiff_t	m,
		std::ptrdiff_t	n,
		double	alpha,
		const double *	a,
		std::ptrdiff_t	lda,
		double *	b,
		std::ptrdiff_t	ldb
	)

inline

B ← α op(A)⁻¹ B (side L) or α B op(A)⁻¹ (R), column-major, B m × n and the triangular A (uplo U or L, its other triangle not read; a diag Unit diagonal not read either) m × m or n × n, as BLAS's dtrsm; A not read where α = 0. Every case is ‘A’ X = B'withA'triangular: Aor a transposed copy of it,Bor a transposed copy of it (X op(A) = Bisop(A)ᵀ Xᵀ = Bᵀ), solved by \ref detail::triangular_solve. A zero on the diagonal gives infinities, as the reference's division does. \return 0, or the position of the first illegal argument, as the reference passes it toXERBLA`.

Definition at line 968 of file blas.hpp.

Variable Documentation

◆ gemv_rows

constexpr std::size_t numerical::blas::gemv_rows {64}

inlineconstexpr

The rows (or, for Aᵀ x, the columns) each task of gemv takes, 64 so that an order of a few hundred shares out (256 left order 300 two tasks).

Definition at line 624 of file blas.hpp.

◆ gemv_weight

constexpr std::size_t numerical::blas::gemv_weight {4}

inlineconstexpr

The weight of gemv's multiply-adds against a compute-bound product's: each reads an element of A once, so the product is bound by the memory, and spreads over more threads' caches for its work (six x86-64 AVX2 cores, Aᵀ x: order 300 17.7 → 7.5 µs, order 1000 74 → 60 µs).

Definition at line 629 of file blas.hpp.

◆ lapack_precision

constexpr double numerical::blas::lapack_precision {0x1p-53}

inlineconstexpr

LAPACK's relative machine precision, ‘dlamch('E’): the unit roundoff2⁻⁵³, half of std::numeric_limits<double>::epsilon(). The default rank cutoff ofdgelsd` is relative to it.

Definition at line 40 of file blas.hpp.

◆ syrk_block

constexpr std::size_t numerical::blas::syrk_block {128}

inlineconstexpr

The columns of syrk's smallest diagonal blocks, formed whole: their wasted half, at most 64 n k multiply-adds, stays small beside the n² k / 2 of the triangle.

Definition at line 44 of file blas.hpp.

◆ syrk_kept_limit

constexpr std::size_t numerical::blas::syrk_kept_limit {(std::size_t {32} << 20U) / sizeof(double)}

inlineconstexpr

The most elements of syrk's operand copy a thread keeps between calls (32 MiB).

Definition at line 47 of file blas.hpp.

◆ trsm_block

constexpr std::size_t numerical::blas::trsm_block {64}

inlineconstexpr

The rows of trsm's blocks: within one, each column solved alone; past it, the matrix product kernel.

Definition at line 913 of file blas.hpp.

Namespaces

Functions

Variables

Function Documentation

◆ complex_dot()

◆ dot() [1/2]

◆ dot() [2/2]

◆ geev()

◆ gelsd()

◆ gemm()

◆ gemv()

◆ geqrf()

◆ gesdd()

◆ gesv()

◆ getrf()

◆ getrs()

◆ orgqr()

◆ potrf()

◆ potrs()

◆ syevd()

◆ syrk()

◆ trsm()

Variable Documentation

◆ gemv_rows

◆ gemv_weight

◆ lapack_precision

◆ syrk_block

◆ syrk_kept_limit

◆ trsm_block