blas.hpp File Reference

The column-major semantics of the BLAS and LAPACK routines NumericAL stands in for (ddot, dgemm, dsyrk, dgemv, dtrsm, dgesv, dgetrf, dgetrs, dpotrf, dpotrs, dsyevd, dgesdd, dgeqrf, dorgqr, dgelsd, dgeev), as their reference implementations define them, over NumericAL's kernels: what the optional numerical_blas library exports with the Fortran calling convention (blas/numerical_blas.cpp), so that a program linked against BLAS and LAPACK, or a switch such as FlexiBLAS, can run on NumericAL. More...

#include <algorithm>
#include <array>
#include <atomic>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <limits>
#include <optional>
#include <utility>
#include <vector>
#include "decomposition.hpp"
#include "eigen.hpp"
#include "gemm.hpp"
#include "matrix.hpp"
#include "parallel.hpp"

Include dependency graph for blas.hpp:

Go to the source code of this file.

Classes
struct	numerical::blas::detail::syrk_operands
	The operands of syrk_triangle, which are row-major P = op(A) and Q = op(A)ᵀ (one of them scaled by α), the column-major C and which triangle of it is written. More...

Namespaces
namespace	numerical
	The NumericAL public API.
namespace	numerical::blas
namespace	numerical::blas::detail

Functions
buffer< double > &	numerical::blas::detail::syrk_operand () noexcept
	The storage of syrk's operand copy this thread keeps between calls.
constexpr bool	numerical::blas::detail::is_operation (char op) noexcept
	Whether op is one of BLAS's N, T or C (either case): no transpose, transpose, conjugate one.
constexpr bool	numerical::blas::detail::transposes (char op) noexcept
	Whether op transposes (T or C: a real matrix's conjugate transpose is its transpose).
constexpr bool	numerical::blas::detail::conjugates (char op) noexcept
	Whether op is BLAS's C: the conjugate transpose (a real matrix's is its transpose).
template<typename T >
bool	numerical::blas::detail::finite (const T &x)
	Whether the real x is finite (an overload per kind of scalar, as conj_of: each one expression).
template<typename U >
bool	numerical::blas::detail::finite (const std::complex< U > &x)
	Whether both parts of the complex x are finite.
template<typename T >
void	numerical::blas::detail::row_major_into (buffer< T > &out, const T *a, std::ptrdiff_t lda, std::ptrdiff_t rows, std::ptrdiff_t cols, T scale, bool conjugate=false)
	The column-major rows × cols a (lda apart) into the row-major out (resized to it), by tiles, each element times scale (conjugated first where conjugate).
template<typename T >
buffer< T >	numerical::blas::detail::row_major (const T *a, std::ptrdiff_t lda, std::ptrdiff_t rows, std::ptrdiff_t cols, T scale, bool conjugate=false)
	rows × cols of the column-major a (lda apart), as a row-major buffer, each element times scale (conjugated first where conjugate).
template<typename T >
buffer< T > &	numerical::blas::detail::kept_storage (std::size_t slot) noexcept
	Row-major storage of the scalar T this thread keeps between calls, in slot (0 for a call's A, 1 for its B): a fresh one of megabytes is mapped anew and faulted in page by page on every call (a solve of order 1000 spent 2.5 ms of its 30 on each).
template<typename T >
buffer< T >	numerical::blas::detail::take_kept (std::size_t slot)
	kept_storage's slot, taken over by the caller (left empty there meanwhile).
template<typename T >
void	numerical::blas::detail::give_back (std::size_t slot, buffer< T > &&storage)
	storage given back to kept_storage's slot, unless past kept_bytes.
template<typename T >
bool	numerical::blas::detail::is_identity (const T *b, std::size_t ldb, std::size_t n)
	Whether the column-major n × n b (ldb apart) is the identity to the bit (a -0 is no zero of it, a solve from it signing zeros of its own): column j against n elements of a unit column of 2n − 1, its one in the middle, from n − 1 − j, a comparison of bytes a column. The first column alone, most other right-hand sides told apart there; the rest by columns on the threads, bound by the memory.
template<typename T >
void	numerical::blas::detail::column_major_into (const T *from, std::size_t rows, std::size_t cols, T *a, std::size_t lda)
	The row-major rows × cols from into the column-major a (lda apart), by tiles.
template<typename T >
void	numerical::blas::detail::column_major_permuted_into (const T *from, std::size_t n, const std::size_t *to, T *a, std::size_t lda)
	column_major_into for the n × n from, its column j written as column to[j] of a.
template<typename T >
void	numerical::blas::detail::product_by_op_a (char transa, std::size_t m, std::size_t n, std::size_t k, const T *l, std::size_t ldl, const T *a, std::size_t lda, T *c, std::size_t ldc)
	C += L op(A)ᵀ (C row-major m × n, L m × k), a the column-major op(A)'s storage (lda apart): op(A)ᵀ is A read in place where transa does not transpose, else a copy of its transpose (an overload per scalar, as conj_of).
void	numerical::blas::detail::product_by_op_a (char transa, std::size_t m, std::size_t n, std::size_t k, const double *l, std::size_t ldl, const double *a, std::size_t lda, double *c, std::size_t ldc)
	product_by_op_a for double: a transposed A read in place too, as op(A)ᵀ's transpose (numerical::detail::gemm_bt, the bits of the product by the copy; NumPy's a @ b.T).
template<typename T >
void	numerical::blas::detail::lower_triangle_rows (T *a, std::size_t lda, std::size_t n, T *rows, bool in)
	The lower triangle of the column-major n × n a (lda apart) into that of the row-major rows (n apart) where in, out of it back into a otherwise, by tiles, bands of them on the threads; the other triangles untouched.
template<typename T >
buffer< T >	numerical::blas::detail::columns_as_rows (const T *a, std::ptrdiff_t lda, std::ptrdiff_t rows, std::ptrdiff_t cols, T scale)
	The cols columns of rows elements of the column-major a (lda apart) as the rows of a row-major buffer, each element times scale: aᵀ, row-major.
template<typename T >
void	numerical::blas::detail::product_of_ops (char transa, char transb, T alpha, std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc)
	C += α op(B)ᵀ op(A)ᵀ, row-major C (m × n, ldc apart), the column-major C of BLAS's gemm as its transpose. A column-major matrix is row-major as its transpose, so op(B)ᵀ is B read in place when B is not transposed (unless α scales a copy), and op(A)ᵀ is A read in place when A is not; a transposed one is copied (an overload per scalar, as conj_of).
void	numerical::blas::detail::product_of_ops (char transa, char transb, double alpha, std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, const double *b, std::size_t ldb, double *c, std::size_t ldc)
	product_of_ops for double: with α = 1 a transposed B is read in place too, as op(B)ᵀ's transpose (numerical::detail::gemm_at, or numerical::detail::gemm_atbt with a transposed A: the bits of the products by the copies; NumPy's a.T @ b).
template<typename Index >
void	numerical::blas::detail::interchanges (const std::vector< std::size_t > &pivots, Index *ipiv)
	The permutation pivots (row i of P A is row pivots[i] of A) as LAPACK's successive row interchanges, 1-based, in ipiv: at step i, row i was interchanged with row ipiv[i].
template<typename T , typename Index >
int	numerical::blas::detail::rectangular_lu (std::size_t m, std::size_t n, T *a, std::size_t lda, Index *ipiv)
	The LU factorisation with partial pivoting of a rectangular column-major A in place, as LAPACK's dgetf2, column by column: NumericAL's blocked LU takes square matrices.
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.
constexpr std::size_t	numerical::blas::detail::strided (std::size_t i, std::size_t length, std::ptrdiff_t inc) noexcept
	Index of element i of a BLAS vector of length elements inc apart: from the last one where inc is negative, as the reference walks it.
template<typename T >
T	numerical::blas::detail::real_dot (std::ptrdiff_t n, const T *x, std::ptrdiff_t incx, const T *y, std::ptrdiff_t incy)
	numerical::blas::dot of the real scalar T.
double	numerical::blas::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	numerical::blas::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 >	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.
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)
	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.
void	numerical::blas::detail::syrk_product (const syrk_operands &operands, std::size_t c0, std::size_t cols, std::size_t r0, std::size_t rows, double *out, std::size_t ldo)
	P(c0…) Q(…, r0…) added to the cols × rows row-major block at out (ldo apart): numerical::detail::gemm, or, where Q is Pᵀ, numerical::detail::gemm_bt from P's rows (the same bits, no copy of Q).
void	numerical::blas::detail::syrk_triangle (syrk_operands operands, std::size_t j0, std::size_t w)
	Adds P Q to the triangle of the diagonal block [j0, j0 + w) of C: halved until it is a block of syrk_block columns or fewer, formed whole aside; above that, the off-diagonal half is one product of (w / 2)² entries, which the matrix product kernel runs on every core. A block of C (rows r0…, columns c0…), column-major, is row-major Cᵀ: it adds P(c0…) Q(…, r0…).
void	numerical::blas::detail::syrk_all (const syrk_operands &operands, std::size_t n)
	The whole triangle of C (order n) from operands: a product the matrix product kernel packs runs in one pass over the triangle's tiles (numerical::detail::gemm_triangle; each operand read where it lies, as syrk_product reads it), a smaller one by halves (syrk_triangle). Column-major C is row-major Cᵀ = P Q, so its upper triangle is the product's lower one.
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)
	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.
void	numerical::blas::detail::triangular_solve (bool lower, bool unit, std::size_t k, std::size_t n, const double *a, std::size_t lda, double *b, std::size_t ldb)
	A X = B in place for the column-major triangular k × k a (lower or upper, a unit diagonal not read) and k × n b: by blocks of trsm_block rows, from the top (lower) or the bottom (upper); a block's diagonal triangle column by column of B, the columns in parallel, then its terms out of the rows still to solve through gemm.
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)
	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`.
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.
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.
template<typename T >
int	numerical::blas::detail::column_cholesky (bool lower, std::size_t size, T *a, std::size_t ld)
	The Cholesky factor of the column-major order-size A in place, column by column as the reference's dpotf2 computes it (lower: L below the diagonal, else U above it, as conjugates for a complex A): where it holds a NaN or an infinity, which NumericAL's Cholesky refuses, the NaN carried through.
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.
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).
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.
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).
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.
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)
	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	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)
	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	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.
void	numerical::blas::detail::real_form (const matrix< std::complex< double > > &v, const vector< std::complex< double > > &values, const std::vector< std::size_t > &column_of, bool conjugate, double *out, std::size_t ld)
	The eigenvectors v (columns, complex) in dgeev's real form, into the column-major out: column j of a real eigenvalue its vector; of a pair λ_j, λ_{j+1} = conj(λ_j) (positive imaginary part first) the real and imaginary parts of λ_j's vector; conjugated where conjugate (the left vectors).
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)
	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	numerical::blas::lapack_precision {0x1p-53}
	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.
constexpr std::size_t	numerical::blas::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	numerical::blas::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	numerical::blas::detail::transpose_tile {32}
	The tile side of the transposing copies between column-major arguments and row-major storage: a tile's rows and columns both stay in the cache, where a whole column read across rows lda apart misses on every element (a 1000 × 1000 copy took several milliseconds).
constexpr std::size_t	numerical::blas::detail::transpose_weight {4}
	The weight of an element of the transposing copies against a multiply-add: bound by the memory, a copy spreads over threads from a few hundred thousand elements (an inverse of order 1000 on six AVX2 cores spent a fifth of its time in its four copies on one thread). The bands capture their bounds by value: by reference, the stores might have changed them, and each was read again (order 100: 5.9 → 10 µs).
constexpr std::size_t	numerical::blas::detail::kept_bytes {std::size_t {32} << 20U}
	The most bytes a thread keeps in each of kept_storage's slots between calls (32 MiB).
constexpr std::size_t	numerical::blas::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	numerical::blas::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	numerical::blas::trsm_block {64}
	The rows of trsm's blocks: within one, each column solved alone; past it, the matrix product kernel.

Detailed Description

The column-major semantics of the BLAS and LAPACK routines NumericAL stands in for (ddot, dgemm, dsyrk, dgemv, dtrsm, dgesv, dgetrf, dgetrs, dpotrf, dpotrs, dsyevd, dgesdd, dgeqrf, dorgqr, dgelsd, dgeev), as their reference implementations define them, over NumericAL's kernels: what the optional numerical_blas library exports with the Fortran calling convention (blas/numerical_blas.cpp), so that a program linked against BLAS and LAPACK, or a switch such as FlexiBLAS, can run on NumericAL.

Each function checks its arguments as the reference does and returns the position of the first illegal one (the number the reference passes to XERBLA), negated for the LAPACK routines as their INFO; it computes nothing then. A matrix holding a NaN or an infinity, which NumericAL's decompositions refuse, is an illegal value of that argument.

Definition in file blas.hpp.