simd.hpp File Reference

The native vector of float and double the kernels hold their lanes in (GNU vector extensions, under GCC and clang; a kernel keeps a scalar form beside it for other compilers). More...

#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstring>
#include <type_traits>

Include dependency graph for simd.hpp:

Go to the source code of this file.

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

Concepts
concept	numerical::detail::split_complex
	A complex number of float or double, which gemm_split multiplies through the real kernel and the kernels take apart into its real and imaginary parts.

Macros
#define	NUMERICAL_ALWAYS_INLINE   [[gnu::always_inline]] inline

Functions
template<typename T >
std::array< T, symmetric_row_group >	numerical::detail::symmetric_rows_tail (const T *rows, std::size_t n, const T *v, T *out, std::size_t i, std::size_t done, std::array< T, symmetric_row_group > sums)
	The columns done…i+r of symmetric_rows, one by one: for each row r in order, out[k] += row_r[k] · v[i+r] for done ≤ k < i + r and row_r[k] · v[k] added to sums[r], its diagonal term last.
template<std::size_t Rows, std::size_t Cols, typename T , typename Full >
void	numerical::detail::through_edge (T *c, std::size_t ldc, std::size_t rows, std::size_t cols, Full full) noexcept
	One tile of a matrix product kernel, rows × cols of C (at most Rows × Cols, leading dimension ldc), through full, which computes a whole tile in place given its first element and leading dimension: a whole tile in C itself, a partial one in a copy whose outside is zero. Only that copy is zeroed: zeroing it for every tile, whole ones included, took a twelfth of an x86-64 product's time.
template<typename T >
std::size_t	numerical::detail::subtract_scaled_vectors (T *, const T *, const T &, std::size_t)
	No leading elements in vectors, for a type no native vector holds.
template<typename T , typename Magnitude >
std::size_t	numerical::detail::largest_magnitude_vectors (const T *, std::size_t, Magnitude &, std::size_t &)
	Only x[0], the one best and at come in with, for a type no native vector holds.
template<typename T >
std::size_t	numerical::detail::subtract_combined_vectors (T *, const T *, const T *, std::size_t, std::size_t, std::size_t)
	No leading blocks in vectors, for a type no native vector holds or a compiler that vectorises the scalar block itself (Clang).
template<typename T >
std::size_t	numerical::detail::subtract_rows_vectors (T *, std::size_t, const T *, std::size_t, const T *, std::size_t, std::size_t, std::size_t)
	No leading columns in vectors, for a type no native vector holds or a compiler that vectorises the scalar block itself (Clang).
template<typename T >
std::size_t	numerical::detail::add_scaled_vectors (T *, const T *, const T &, std::size_t)
	No leading elements in vectors, for a type no native vector holds.
template<typename T >
T	numerical::detail::symmetric_row (const T *row, const T *v, const T &f, T *out, std::size_t count)
	Σ row[i] · v[i] and out[i] += row[i] · f, in order, for a type no native vector holds.
template<typename T >
std::array< T, symmetric_row_group >	numerical::detail::symmetric_rows (const T *rows, std::size_t n, const T *v, T *out, std::size_t i)
	symmetric_rows one column at a time, for a type no native vector holds.
template<typename T >
T	numerical::detail::dot_product (const T *x, const T *y, std::size_t count)
	Σ x[i] · y[i], in order, for a type no native vector holds.
template<typename T >
T	numerical::detail::conjugate_dot_product (const T *x, const T *y, std::size_t count)
	Σ conj(x[i]) · y[i], in order, where no native vector holds the parts.
template<typename T >
void	numerical::detail::row_dots (const T *m, std::size_t ld, std::size_t rows, const T *x, std::size_t count, T *out)
	dot_product of each row of m with x, for a type no native vector holds.
template<typename T >
T	numerical::detail::operator_times (const T &a, const T &b)
	a · b, as std::operator* computes it, in line. libc++'s complex product (the four products and their sums as statements of their own, so unfused) is copied here, its rescue of a product whose parts are both NaN left to it: where it was a call, out of line in a translation unit as large as the Python extension, it took 45 % of a complex LU of order 48. Elsewhere (libstdc++, whose product the compiler expands) and for a real T, the operator itself.
template<typename T >
void	numerical::detail::subtract_product (T *y, const T &f, const T *x, std::size_t count)
	y[i] -= f · x[i] for count elements of y and x, which do not overlap. Each element is computed exactly as that scalar statement computes it (its product rounded before its difference, or fused with it where the compiler contracts the statement); the arguments follow the expression, since a compiler may contract f · x and x · f of complex numbers differently. Under GCC, over native vectors where T has them: at -O2 GCC leaves the scalar loop scalar, unable to rule out an overlap without a check it will not add. Clang vectorises the scalar loop itself, with that check, and better than these vectors (an LU solve of order 128 on arm64 took 19 % longer through them). For double, gemm.hpp's overload chooses the x86-64 units at run time, a fused multiply-add per element there.
template<typename T >
void	numerical::detail::subtract_products (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count)
	y[i] -= Σₜ f[t] · x[(t · ldx) + i] for count elements of y, terms rows of x apart by ldx: subtract_product for each term in ascending order, so each element is computed as that sequence of calls computes it. For double, gemm.hpp's overload holds the elements in registers through the terms.
template<typename T > requires std::is_floating_point_v<T>
void	numerical::detail::subtract_products (T *y, const T *f, const T *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_products for a real T: eight elements of y held in a local array through the terms, which the compiler keeps in registers, where each term was a call through the whole row. Each element takes its terms in the same order as a sequence of subtract_product gives them (an LU solve of order 64 and four columns on arm64: 18 → 8.5 µs), and the same value under Clang, and under GCC on x86-64, which fuses neither; under GCC on arm64 the vectors of subtract_combined_vectors and subtract_scaled_vectors round apart and the scalar statements fused, and their splits differ.
template<typename T > requires std::is_floating_point_v<T>
void	numerical::detail::subtract_rows (T *y, std::size_t ldy, const T *f, std::size_t ldf, const T *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_products for product_rows rows of y at once, ldy apart, each with its own factors (ldf apart in f) on the same rows of x: y_r[i] -= Σₜ f_r[t] · x[(t · ldx) + i]. Each element takes the statement of subtract_products, a term at a time in ascending order, so the same value; each element of x read serves every row, and the rows' chains of multiply-adds advance together, where a row alone waits on its last one through each term.
template<typename T >
void	numerical::detail::subtract_product (T *y, const T *x, const T &f, std::size_t count)
	y[i] -= x[i] · f, as the other subtract_product (a real product is the same either way round).
template<split_complex C>
void	numerical::detail::complex_scaled_add (C &y, const C &f, const C &x)
	y += f · x of one complex number, its part products rounded apart as libc++'s operator computes them (add_scaled_vectors's arithmetic, for the elements past its vectors).
template<split_complex C>
void	numerical::detail::add_product (C *y, const C &f, const C *x, std::size_t count)
	add_product of complex numbers: vectors of their parts, then complex_scaled_add, every element the same arithmetic.
template<typename T >
void	numerical::detail::add_one_product (T &y, const T &f, const T &x)
	y += f · x of one element, as add_product computes each of its elements.
template<split_complex C>
void	numerical::detail::add_one_product (C &y, const C &f, const C &x)
	add_one_product of complex numbers: complex_scaled_add, as add_product computes them.
template<typename T >
void	numerical::detail::add_product (T *y, const T &f, const T *x, std::size_t count)
	y[i] += f · x[i], as subtract_product.

Variables
constexpr std::size_t	numerical::detail::transpose_tile {32}
	The side of the tiles a transposing copy goes by (an einsum permutation of two axes, the eigenvectors' columns): a tile read and a tile written fit the L1 cache, where element by element each write, a row apart, missed it (order 1024 on arm64: 4.7 → 1.5 ms; the eigenvectors of order 1024, a tenth of an eigh on six x86-64 cores).
constexpr std::size_t	numerical::detail::symmetric_row_group {4}
	The rows of a lower triangle symmetric_rows takes at once.
constexpr std::size_t	numerical::detail::product_rows {8}
	The rows of y subtract_rows takes at once.

Detailed Description

The native vector of float and double the kernels hold their lanes in (GNU vector extensions, under GCC and clang; a kernel keeps a scalar form beside it for other compilers).

Definition in file simd.hpp.

Macro Definition Documentation

NUMERICAL_ALWAYS_INLINE

#define NUMERICAL_ALWAYS_INLINE   [[gnu::always_inline]] inline

Definition at line 22 of file simd.hpp.