gemm.hpp File Reference

The matrix product kernel behind matrix * matrix and the einsum matmul path. More...

#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <atomic>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
#include "parallel.hpp"
#include "amx.hpp"
#include "buffer.hpp"
#include "simd.hpp"
#include "x86.hpp"

Include dependency graph for gemm.hpp:

Go to the source code of this file.

Classes
struct	numerical::detail::b_rows
	B as gemm_packed reads it by default: row-major, k × n, its rows ldb apart. More...
struct	numerical::detail::b_transposed
	B given as Bᵀ: row-major, n × k, its rows ldb apart (C = A·Bᵀᵀ, no copy of the transpose). More...
struct	numerical::detail::a_rows
	A as gemm_packed reads it by default: row-major, m × k, its rows lda apart, packed by the tile. More...
struct	numerical::detail::a_transposed
	A given as Aᵀ: row-major, k × m, its rows lda apart (C = Aᵀᵀ·B, no copy of the transpose). More...
struct	numerical::detail::vector_tile_kernel< T >
	The vector path's tile: gemm_mr × gemm_nr elements of C in registers (gemm_kernel), each product rounded before its sum unless the compiler contracts c + a·b into a fused multiply-add. More...
struct	numerical::detail::vector_tile_kernel< T >::session
	Nothing to turn on around a task. More...
class	numerical::detail::packing_storage< T >
	The storage of a product's packed operands, kept by the calling thread from one product to the next. More...
struct	numerical::detail::split_operands< R >
	The real operands of gemm_split, kept by each thread from one product to the next (split_storage): mapping their pages anew took a third to a half of a complex product of order 64 to 256 on one x86-64 core (six AVX2 cores, all threads: order 128 0.52 → 0.17 ms), 15–20 % of one of order 256 on arm64. More...
struct	numerical::detail::packed_left< T >
	A left operand packed once for several products by it (an LU's −L₂₁ against blocks of its U₁₂): the tile chosen and A packed for it (gemm_packed's layout of one panel), or A to read in place (plain). More...

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

Macros
#define	NUMERICAL_OUT_OF_LINE   [[gnu::noinline]]
#define	NUMERICAL_SCALAR_DEPTH
#define	NUMERICAL_UNROLL_BLOCK

Typedefs
template<typename T >
using	numerical::detail::gemm_tile = std::array< T, gemm_mr *gemm_nr< T > >
template<term_rounding Mode>
using	numerical::detail::rounding_is = std::integral_constant< term_rounding, Mode >
	The rounding mode Mode as a type, which picks rounded_term's overload.

Enumerations
enum class	numerical::gemm_path : std::uint8_t {   numerical::portable = 0 , numerical::neon = 1 , numerical::avx2 = 2 , numerical::avx512 = 3 ,   numerical::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...
enum class	numerical::detail::term_rounding : unsigned char { numerical::detail::compiler , numerical::detail::fused , numerical::detail::x86_fused }
	How numerical::detail::gemm_rows rounds each term sum + a·b, as the tile of the path in use does: as the compiler rounds it (the vector tile's way), one fused multiply-add (arm64's NEON and AMX tiles), or one fused multiply-add of AVX2's unit (the x86-64 tiles), in code built for any x86-64 processor. More...
enum class	numerical::detail::left_kernel : unsigned char {   numerical::detail::plain , numerical::detail::vector , numerical::detail::neon , numerical::detail::amx ,   numerical::detail::avx512 , numerical::detail::avx2 }
	The tile a product's left operand is packed for (numerical::detail::pack_left): plain where it is not packed. More...

Functions
bool	numerical::detail::gemm_path_usable (gemm_path path) noexcept
	Whether path can run here.
gemm_path	numerical::detail::best_gemm_path (gemm_path most) noexcept
	The most capable path that can run here, most at the highest.
std::atomic< gemm_path > &	numerical::detail::gemm_path_cap () noexcept
	The most capable path products may take, as numerical::limit_gemm_path set it: none is probed until a product (or a question) reaches it, so a cap below the AMX keeps its instructions from ever running.
std::atomic< std::uint8_t > &	numerical::detail::gemm_path_found () noexcept
	The path gemm_path_current found last, packed with the cap it was found under (cap · 8 + path); 0xFF before any.
gemm_path	numerical::detail::gemm_path_current () noexcept
	The path the next double product takes: the most capable one under the cap that runs here (each path's probe is made once). Kept with its cap, one load, where each call walked the paths down from the cap through their probes' guards: every small step of a decomposition asks (an LU of order 100 through dgetrf_ on one AVX2 core spent 5 % of its time here). A race stores the same value twice.
void	numerical::detail::subtract_product (double *y, const double &f, const double *x, std::size_t count)
	y[i] −= f · x[i] for double on the widest x86-64 units the matrix product's kernel may use (gemm_path_current), each element one fused multiply-add there; elsewhere subtract_product's template. The decompositions' row updates take it: a wheel built for the baseline x86-64 holds two double per vector and fuses none.
void	numerical::detail::subtract_product (double *y, const double *x, const double &f, std::size_t count)
	y[i] −= x[i] · f for double, as the other (a real product is the same either way round).
void	numerical::detail::subtract_products (double *y, const double *f, const double *x, std::size_t ldx, std::size_t terms, std::size_t count)
	y[i] -= Σₜ f[t] · x[(t · ldx) + i] for double, the units chosen as subtract_product chooses them, the elements held in registers through the terms there.
void	numerical::detail::subtract_rows (double *y, std::size_t ldy, const double *f, std::size_t ldf, const double *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_rows for double, the units chosen as subtract_products chooses them: each element the value subtract_products gives it.
void	numerical::detail::subtract_products (std::complex< double > *y, const std::complex< double > *f, const std::complex< double > *x, std::size_t ldx, std::size_t terms, std::size_t count)
	subtract_products for std::complex<double>: on AVX-512 or AVX2 with FMA where the matrix product's kernel runs there, each part of each term one fused multiply-add (Re f · x, then Im f times the swapped parts), the elements in registers through the terms; elsewhere the template, whole complex products as the compiler rounds them.
void	numerical::detail::subtract_product (std::complex< double > *y, const std::complex< double > &f, const std::complex< double > *x, std::size_t count)
	y[i] −= f · x[i] for std::complex<double>: one term of subtract_products on AVX-512 or AVX2 where the matrix product's kernel runs there, else the template.
void	numerical::detail::subtract_product (std::complex< double > *y, const std::complex< double > *x, const std::complex< double > &f, std::size_t count)
	y[i] −= x[i] · f for std::complex<double>, as the other (the template keeps the expression's order).
void	numerical::detail::add_product (double *y, const double &f, const double *x, std::size_t count)
	y[i] += f · x[i] for double, as subtract_product.
void	numerical::detail::add_product (std::complex< float > *y, const std::complex< float > &f, const std::complex< float > *x, std::size_t count)
	y[i] += f · x[i] for complex float: on the x86-64 paths by AVX2 vectors of its arithmetic (avx2_scaled_add_c32), the same bits as add_product's template, which runs elsewhere.
void	numerical::detail::add_product (std::complex< double > *y, const std::complex< double > &f, const std::complex< double > *x, std::size_t count)
	add_product for complex double (avx2_scaled_add_c64 on the x86-64 paths).
double	numerical::detail::dot_product (const double *x, const double *y, std::size_t count)
	Σ x[i] · y[i] for double, as subtract_product chooses its units (in another order of summation there, as fixed).
template<split_complex C>
C	numerical::detail::complex_dot_product (const C *x, const C *y, std::size_t count)
	Σ x[i] · y[i] of complex numbers: on the x86-64 units with AVX2 (the AVX-512 path too) their fused chains (avx2_complex_dot), elsewhere the native vectors' (in another order of summation there, as fixed).
std::complex< double >	numerical::detail::dot_product (const std::complex< double > *x, const std::complex< double > *y, std::size_t count)
	dot_product for std::complex<double>, through complex_dot_product.
std::complex< float >	numerical::detail::dot_product (const std::complex< float > *x, const std::complex< float > *y, std::size_t count)
	dot_product for std::complex<float>, through complex_dot_product.
template<split_complex C>
void	numerical::detail::complex_row_dots (const C *m, std::size_t ld, std::size_t rows, const C *x, std::size_t count, C *out)
	row_dots of complex numbers: each row by complex_dot_product, as chunked_dot sums a row.
void	numerical::detail::row_dots (const std::complex< double > *m, std::size_t ld, std::size_t rows, const std::complex< double > *x, std::size_t count, std::complex< double > *out)
	row_dots for std::complex<double>, through complex_row_dots.
void	numerical::detail::row_dots (const std::complex< float > *m, std::size_t ld, std::size_t rows, const std::complex< float > *x, std::size_t count, std::complex< float > *out)
	row_dots for std::complex<float>, through complex_row_dots.
void	numerical::detail::row_dots (const double *m, std::size_t ld, std::size_t rows, const double *x, std::size_t count, double *out)
	row_dots for double: on the x86-64 units each row by their dot_product (wider than four rows of the portable vectors), elsewhere four rows at a time.
template<typename T >
void	numerical::detail::gemm_kernel (std::size_t kc, const T *ap, const T *bp, std::size_t ldb, gemm_tile< T > &tile)
	tile += Aₚ·Bₚ over one panel, for any scalar.
template<std::size_t Mr, typename T >
void	numerical::detail::gemm_pack_a (const T *a, std::size_t lda, std::size_t rows, std::size_t kc, T *out)
	Copies rows × kc of A (from row-major, leading dimension lda) into slivers of Mr rows, each stored column by column; a last partial sliver is padded with zeros.
template<std::size_t Nr, typename T >
void	numerical::detail::gemm_pack_b (const T *b, std::size_t ldb, std::size_t kc, std::size_t cols, T *out)
	Copies kc × cols of B (from row-major, leading dimension ldb) into slivers of Nr columns, each stored row by row; a last partial sliver is padded with zeros.
template<std::size_t Nr, typename T >
void	numerical::detail::gemm_pack_bt (const T *bt, std::size_t ldbt, std::size_t kc, std::size_t cols, T *out)
	gemm_pack_b of B = Bᵀᵀ, from the row-major bt (cols × kc, leading dimension ldbt): each column of a sliver is a row of bt, read along it; a last partial sliver is padded with zeros.
template<typename T >
const T *	numerical::detail::rows_in_place (const T *x, std::size_t ld, std::size_t &ld_out) noexcept
	A row-major rows × cols operand for the plain loop, x itself (ld apart).
template<typename T >
const T *	numerical::detail::rows_from_transpose (const T *x, std::size_t ld, std::size_t rows, std::size_t cols, buffer< T > &scratch, std::size_t &ld_out)
	A row-major rows × cols operand for the plain loop, copied into scratch from its transpose x (cols × rows, ld apart).
template<std::size_t Mr, typename T >
void	numerical::detail::gemm_pack_at (const T *at, std::size_t ldat, std::size_t rows, std::size_t kc, T *out)
	The slivers of Mr rows of A (rows × kc), column by column, from its transpose at (kc × rows row-major, ldat apart): each column of a sliver is a run of a row of at; a last partial sliver is padded with zeros, as gemm_pack_a lays them out.
template<typename T >
T	numerical::detail::rounded_apart_add (T sum, T a, T b) noexcept
	sum + a·b as the vector tile rounds it (vector_tile_kernel), so that the plain loop's elements are its own: under Clang's -ffp-contract=on one expression, fused in both; under GCC's, whose vector extensions fuse nothing while a scalar statement's multiply-add fuses, the product apart, a statement of its own.
template<typename R >
std::complex< R >	numerical::detail::rounded_apart_add (std::complex< R > sum, std::complex< R > a, std::complex< R > b) noexcept
	rounded_apart_add of complex numbers: the product in line (operator_times, the operator's own rounding, unfused), where it was a call that took a sixth of a complex LU.
template<typename T >
T	numerical::detail::rounded_term (rounding_is< term_rounding::compiler >, T sum, T a, T b) noexcept
	sum + a·b as the compiler rounds it (rounded_apart_add).
template<typename T >
T	numerical::detail::rounded_term (rounding_is< term_rounding::fused >, T sum, T a, T b) noexcept
	sum + a·b as one fused multiply-add.
template<std::size_t R, std::size_t N, term_rounding Mode, typename T >
void	numerical::detail::narrow_block (std::size_t k, const T *a, std::size_t lda, const T *b, std::size_t ldb, T *c, std::size_t ldc) noexcept
	C += A·B for R rows and N columns of C, each sum in ascending order of the depth, as the plain loop's: the block held in registers, of a size the compiler knows; each term rounded as Mode says. Out of line: its unrolled sizes, inlined, spent GCC's growth budget, and the LU factorisation, which runs none of them, lost 3 to 8 % at every order (one arm64 thread, GCC 14).
template<std::size_t R, term_rounding Mode, typename T >
void	numerical::detail::narrow_rows (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) noexcept
	C += A·B for R whole rows of C of n columns: blocks of narrow_columns, the last narrower.
template<term_rounding Mode, typename T >
void	numerical::detail::gemm_rows_rounded (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 += A·B by the plain i-k-j loop, parallel over the rows of C. Past a depth of 1 (narrow_least_depth), the rows run narrow_group at a time by blocks (narrow_rows), the last rows one by one: the same sums, six to ten times as fast from order 4 to 16 (M1 Pro) as the loop of unknown width, whose sums the compiler cannot keep in registers. A depth of 1 streams as fast through that loop.
bool	numerical::detail::rows_fuse () noexcept
	Whether the tile a double or float product past the plain loop runs on fuses each term (gemm_rows rounds its own so): double on arm64's NEON and AMX paths, float on the NEON tile there. The x86-64 tiles' double takes term_rounding::x86_fused; the vector tile, and so the portable path, rounds as the compiler does, as the plain loop then does. On arm64 every cap but the portable one runs the NEON or the AMX path, so the cap alone answers: one load, where the path in use took an order-4 product on one arm64 core from 10.6 to 13.3 ns (the cap, 11.5).
template<typename T >
void	numerical::detail::gemm_rows (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 += A·B by the plain i-k-j loop (gemm_rows_rounded), each element rounded as the tile of the path in use rounds it, so that it is the same whatever the shape of the product it falls in: a caller that splits its columns by the threads (an LU inverse) then has the same bits on any number of them.
template<typename T >
T *	numerical::detail::line_aligned (T *at) noexcept
	The first element of at's run that starts a cache line (or at's nearest element towards one, for an element whose size does not divide the allocator's alignment).
template<typename Kernel >
constexpr bool	numerical::detail::gemm_shares (std::size_t m, std::size_t n, std::size_t k) noexcept
	Whether the product's threads pay: below Kernel's serial_below multiply-adds, the product runs on the caller alone.
template<typename Kernel , typename T >
void	numerical::detail::prefetch_next_row (const T *, std::size_t) noexcept
	Nothing, for a tile that prefetches no C (prefetches_next_row_of_c).
template<typename Kernel , typename T > requires prefetches_next_row_of_c<Kernel>
void	numerical::detail::prefetch_next_row (const T *c, std::size_t ldc) noexcept
	The C tile at c prefetched by Kernel's prefetch_rows.
template<typename Kernel >
constexpr std::size_t	numerical::detail::task_rows_cap (std::size_t) noexcept
	The most rows of a task of Kernel over a panel of depth kc: its block of rows.
template<typename Kernel > requires requires(std::size_t kc) { Kernel::task_rows_for(kc); }
constexpr std::size_t	numerical::detail::task_rows_cap (std::size_t kc) noexcept
	task_rows_cap for a tile whose task rows depend on the panel's depth (Kernel::task_rows_for).
template<typename Kernel , typename T , typename B = b_rows, typename A = a_rows>
void	numerical::detail::gemm_packed (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, const T *a_packed=nullptr)
	The packed product with the tile Kernel computes: for each block of B's columns and each panel of the depth, both operands are packed (in parallel, by the blocks the product takes), then tasks of gemm_mc rows by gemm_slivers_per_task slivers of columns run in parallel, each writing a disjoint block of C, one tile at a time. Given a_packed (A as pack_left packs it for Kernel, a depth of one panel), A is read from there and only B is packed: the same tiles, the same bits. B is B's layout (b_rows, b_transposed), A A's (a_rows, a_transposed): the packing alone differs, so the bits do not.
template<typename Kernel , typename T , typename B = b_rows, typename A = a_rows>
void	numerical::detail::gemm_packed_triangle (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, bool lower)
	gemm_packed restricted to a triangle of the square C = A·B (n × n): its lower one (j ≤ i) where lower, else its upper one (j ≥ i); the other is neither written nor depended on. The same packing (A's row blocks with no tile in the triangle left out) and the same tiles, so each element's sum is the product's; a tile across the diagonal runs on a copy of its block of C, the triangle's part copied back. A syrk takes the triangle in one pass, its operands packed once a panel, where products by halves packed them again for every block.
template<typename T >
bool	numerical::detail::gemm_elsewhere (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)
	A product a faster path than the vector one takes, if any: none for this scalar.
template<typename B , typename A = a_rows>
bool	numerical::detail::gemm_double_elsewhere (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)
	A double product on the path in use (numerical::gemm_path_in_use), unless it is the portable one, B and A laid out as B and A say (b_rows, a_rows and their transposes).
template<typename A , typename B >
void	numerical::detail::gemm_triangle (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, bool lower)
	The triangle of the double product C = A·B (gemm_packed_triangle) on the tile of the path in use, A and B laid out as A and B say.
bool	numerical::detail::gemm_elsewhere (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)
	gemm_double_elsewhere of a row-major B.
template<typename B , typename A = a_rows>
bool	numerical::detail::gemm_float_elsewhere (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc)
	A float product on the AMX where it suits the product, on the NEON tile on arm64's fused paths, on the x86-64 tiles there; elsewhere the vector tile. B and A laid out as B and A say.
bool	numerical::detail::gemm_elsewhere (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc)
	gemm_float_elsewhere of row-major operands.
template<typename T >
void	numerical::detail::gemm (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 += A·B, every element of C summed in ascending order of the depth.
bool	numerical::detail::split_by_parts (std::size_t m) noexcept
	Whether gemm_split runs a product of m rows of A by four real products of each block's parts (split_parts_block), else by one of its expanded operands (split_expanded_block): the same bits, the parts' products faster where they read B's parts in place (gemm_packed). On the AMX path up to split_parts_rows; elsewhere (the NEON path) up to one row of the NEON tile, its b_in_place_rows: M1 Pro, all threads, by 400 × 500, 1 to 4 rows 2 to 18 % faster, 6 to 32 rows 5 to 52 % slower.
template<typename R >
split_operands< R > &	numerical::detail::split_storage ()
	The calling thread's split_operands. A product uses them whole on one thread, and the real product it runs on them nests no complex one.
template<split_complex Z, typename B = b_rows, typename A = a_rows>
void	numerical::detail::split_parts_block (std::size_t m, std::size_t nb, std::size_t kb, std::size_t pc, std::size_t jc, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, typename Z::value_type *ar, typename Z::value_type *br, typename Z::value_type *cr)
	One block of gemm_split ([Re C | Im C] at cr, m × nb each, += A·B, A m × kb and B kb × nb) by four real products of its parts: Re C += Re A·Re B, then −Im A·Im B; Im C += Re A·Im B, then Im A·Re B, each continuing its elements' chains ((−x)·(−y) is x·y exactly: the bits of −Im A·−Re B). B's parts are a copy of its size, which a product of a few rows of A reads in place. The block is A's columns and B's rows from pc, B's columns from jc, A and B laid out as A and B say.
template<split_complex Z, typename B = b_rows, typename A = a_rows>
void	numerical::detail::split_expanded_block (std::size_t m, std::size_t nb, std::size_t kb, std::size_t pc, std::size_t jc, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, typename Z::value_type *ar, typename Z::value_type *br, typename Z::value_type *cr)
	One block of gemm_split (split_parts_block's) by one real product of its expanded operands, [Re C | Im C] += [Re A | −Im A] · [[Re B, Im B], [Im B, −Re B]]: the same chains, B copied four times its size, the product as wide and as deep as the block's two.
template<split_complex Z, typename B = b_rows, typename A = a_rows>
void	numerical::detail::gemm_split (std::size_t m, std::size_t n, std::size_t k, const Z *a, std::size_t lda, const Z *b, std::size_t ldb, Z *c, std::size_t ldc)
	C += A·B for complex numbers through the real kernel, by blocks of split_columns columns and split_depth of depth, on the real path's units and threads (a scalar complex tile has none). Within a block, each element of C takes the products of the real parts of A in ascending order of the depth, then those of the imaginary parts: by four real products of the block's parts (split_parts_block) or one of its expanded operands (split_expanded_block), as split_by_parts chooses, the same bits either way. A and B laid out as A and B say, read where they lie as the parts are copied.
bool	numerical::detail::plain_loop_pays (std::size_t m, std::size_t n, std::size_t k) noexcept
	Whether a double product of m rows, n columns and depth k, past gemm_rows's threshold, still runs faster in the plain loop than packed on the path in use: packing B's slivers for so few rows of A, or so short a depth, did not pay. One core, the plain loop's time over the packed product's: the NEON tile at a depth of 16 or less whatever the shape (0.5 to 0.6 up to 2000 × 2000), at 64 or less for 64 rows or fewer (0.4 to 0.9; the same under GCC 14 and on eight threads); the AVX2 tile for 16 rows or fewer to a depth of 64, or 32 rows to 32 (0.2 to 0.8, six cores too), as wide as plain_loop_least_width. The AMX's packing pays from 16 rows.
template<typename T >
bool	numerical::detail::plain_loop_suits (std::size_t, std::size_t, std::size_t) noexcept
	plain_loop_pays for a product of T: double's alone, the one it was measured for.
template<>
bool	numerical::detail::plain_loop_suits< double > (std::size_t m, std::size_t n, std::size_t k) noexcept
	plain_loop_pays for double.
template<typename T >
void	numerical::detail::gemm_large (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)
	A product past the plain loop: on the AMX or x86-64 tiles where gemm_elsewhere takes it, else on the vector tile.
template<typename B , typename A , split_complex T>
void	numerical::detail::gemm_large_complex (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)
	A complex product past the plain loop: on the vector tile on the portable path, whole complex products summed in order; on the x86-64 complex tiles there; else (arm64) through the real kernel (gemm_split). B and A laid out as B and A say.
template<typename B , typename A >
void	numerical::detail::gemm_large_layouts (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)
	A double product past the plain loop, B and A laid out as B and A say (gemm_layouts).
template<typename B , typename A >
void	numerical::detail::gemm_large_layouts (std::size_t m, std::size_t n, std::size_t k, const float *a, std::size_t lda, const float *b, std::size_t ldb, float *c, std::size_t ldc)
	A float product past the plain loop, B and A laid out as B and A say (gemm_layouts).
template<typename B , typename A , split_complex T>
void	numerical::detail::gemm_large_layouts (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)
	A complex product past the plain loop, B and A laid out as B and A say (gemm_layouts).
template<typename A , typename B , typename T >
void	numerical::detail::gemm_layouts (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)
	gemm with either operand given as its transpose (A: a_rows or a_transposed, B: b_rows or b_transposed), so that a product by a transpose (a syrk's P Pᵀ, BLAS's transa and transb, a NumPy array in Fortran order) needs no copy of it. The tiles are gemm's, the operands packed (or, through gemm_split, copied into their parts) from where they lie, and so are the bits; a product gemm runs in the plain loop copies a transposed operand row-major and runs there, small as it is.
template<typename T >
void	numerical::detail::gemm_bt (std::size_t m, std::size_t n, std::size_t k, const T *a, std::size_t lda, const T *bt, std::size_t ldbt, T *c, std::size_t ldc)
	gemm_layouts with B given as Bᵀ (bt, n × k, ldbt apart).
template<typename T >
void	numerical::detail::gemm_at (std::size_t m, std::size_t n, std::size_t k, const T *at, std::size_t ldat, const T *b, std::size_t ldb, T *c, std::size_t ldc)
	gemm_layouts with A given as Aᵀ (at, k × m, ldat apart).
template<typename T >
void	numerical::detail::gemm_atbt (std::size_t m, std::size_t n, std::size_t k, const T *at, std::size_t ldat, const T *bt, std::size_t ldbt, T *c, std::size_t ldc)
	gemm_layouts with both operands given as their transposes.
template<typename T >
packed_left< T >	numerical::detail::read_in_place (std::size_t m, std::size_t k, const T *a, std::size_t lda) noexcept
	The m × k a read in place by gemm_with, unpacked: each product packs it as gemm does.
template<typename T >
packed_left< T >	numerical::detail::left_rows (const packed_left< T > &left, std::size_t r0, std::size_t rows) noexcept
	The rows [r0, r0 + rows) of left, r0 a multiple of left.granule: their slivers in its packing, or their rows of A in place.
template<typename T >
left_kernel	numerical::detail::left_kernel_for (std::size_t, std::size_t, std::size_t) noexcept
	The tile gemm runs a product of m × n by depth k on, for a left operand packed once for many products: plain for a scalar other than double (written as an overload per scalar, as conj_of).
template<>
left_kernel	numerical::detail::left_kernel_for< double > (std::size_t m, std::size_t n, std::size_t k) noexcept
	left_kernel_for a double product: past gemm_rows, as gemm_elsewhere chooses, of a depth of one panel; plain otherwise.
template<typename Kernel , typename T >
const T *	numerical::detail::pack_left_for (std::size_t m, std::size_t k, const T *a, std::size_t lda, buffer< T > &storage)
	A packed for Kernel into storage, in parallel by its blocks of rows, as gemm_packed packs it.
template<typename T >
packed_left< T >	numerical::detail::pack_left (std::size_t m, std::size_t, std::size_t k, const T *a, std::size_t lda, buffer< T > &)
	The m × k a packed once for products with it of about n columns each (gemm_with), in storage (grown as needed, kept by the caller from one call to the next): on the tile gemm would choose for m × n × k. A scalar other than double (here), a depth of more than one panel, or a product gemm_rows takes, leaves A to be read in place.
template<>
packed_left< double >	numerical::detail::pack_left< double > (std::size_t m, std::size_t n, std::size_t k, const double *a, std::size_t lda, buffer< double > &storage)
	pack_left for double, on the tile left_kernel_for chooses.
template<typename T >
void	numerical::detail::gemm_with (const packed_left< T > &left, std::size_t n, const T *b, std::size_t ldb, T *c, std::size_t ldc)
	C += A·B for the packed left (pack_left) and the k × n b, as gemm computes it, A read in place here, for a scalar other than double.
template<>
void	numerical::detail::gemm_with< double > (const packed_left< double > &left, std::size_t n, const double *b, std::size_t ldb, double *c, std::size_t ldc)
	gemm_with for double: on the tile A was packed for, A read from its packing, where gemm would choose that tile for this n; else as gemm. Either way gemm's tiles and order of summation, its bits.
template<typename T >
T	numerical::detail::chunked_dot (const T *a, const T *b, std::size_t count)
	Σ a[i] · b[i] for a plain dot product: chunks of dot_chunk elements in parallel, each by dot_product (several vector accumulators, on the widest x86-64 units for double), their sums added in the chunks' order. The same on any number of threads; in the last bits, not the ascending sum, which waits on each addition before the next (ten times slower at 100 000 elements).
template<typename T >
void	numerical::detail::matrix_times_vector (const T *a, std::size_t rows, std::size_t cols, const T *x, T *y)
	y = A x for the row-major rows × cols a: each element the dot product of its row with x, as chunked_dot sums it (a row of a chunk or less is one dot_product), the rows in parallel by groups (row_dots, which share each load of x). The same on any number of threads. A double matrix of amx_dot_least elements to fewer than amx_dot_most, on the AMX path, runs on the AMX (amx_dot_rows_f64): each element then its eight lanes of fused multiply-adds, added in order.
gemm_path	numerical::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	numerical::gemm_path_in_use () noexcept
	The path double matrix products take.

Variables
constexpr bool	numerical::detail::neon_compiled {false}
	Whether the NEON path's tile is compiled in: on every arm64 target GCC or Clang builds.
constexpr std::size_t	numerical::detail::gemm_mr {4}
	The rows of C one register tile covers.
template<typename T >
constexpr std::size_t	numerical::detail::gemm_nr {sizeof(T) >= 16 ? std::size_t {6} : std::size_t {96} / sizeof(T)}
	The columns of C one register tile covers: 96 bytes of elements (12 double, 24 float), which with gemm_mr rows fills the vector registers of both arm64 and x86-64 with AVX.
constexpr std::size_t	numerical::detail::gemm_kc {256}
	The depth of one panel: a B sliver (gemm_kc × gemm_nr) stays in the L1 cache.
constexpr std::size_t	numerical::detail::gemm_mc {64}
	The rows of one task: an A block (gemm_mc × gemm_kc) stays in the L2 cache.
constexpr std::size_t	numerical::detail::gemm_slivers_per_task {4}
	The B slivers of one task.
constexpr std::size_t	numerical::detail::tasks_per_thread {4}
	The tasks a product's compute pass hands each thread at least, while its rows allow: a slower core then holds back no one.
constexpr std::size_t	numerical::detail::gemm_kernel_speedup {16}
	The multiply-adds the register tile does in the time of one in work_per_thread's dependent chain (about 20 G/s on one arm64 or x86-64 AVX2 core): a task weighs its work by it, so a product of 0.1 ms does not start eight threads.
template<typename T >
constexpr std::size_t	numerical::detail::gemm_nc {gemm_nr<T> * 128}
	The columns of B packed at once, a whole number of slivers.
constexpr std::size_t	numerical::detail::gemm_packed_threshold {std::size_t {1} << 13U}
	Below this many multiply-adds the plain loop (gemm_rows, by blocks of rows held in registers) runs instead of the packed kernel, both summing in the same order. Against the packed kernel on 400 random shapes from 512 multiply-adds, the plain loop took 0.25 to 0.41 times as long below 2^13 on arm64 (M1 Pro, geometric means by octave, the slowest shape 0.91); above, it lost up to 2.5 times on shapes the AMX suits. On x86-64 the plain loop has the baseline's vectors and the packed one the AVX2 or AVX-512 units: there it took 0.43 to 0.55 times as long below 2^11, 0.86 above with shapes up to 2 times slower (six AVX2 cores, GCC 13).
constexpr std::size_t	numerical::detail::gemm_min_depth {8}
	Below this depth the plain loop runs too: each panel reads and writes its tiles of C once, which a shallow product does not amortise (1000×1000×1 ran six times slower packed, ×4 still lost 10 %).
constexpr std::size_t	numerical::detail::narrow_least_depth {2}
	The least depth at which rows of C wider than one block (narrow_columns) run by blocks in gemm_rows; one block or narrower takes any depth past 1. Against the plain loop of unknown width, on 500 random shapes below the packed kernel, the blocks took 0.11 to 0.63 times as long from a depth of 2 on arm64 (M1 Pro, geometric means by width); on x86-64, where the loop has the baseline's vectors, 0.35 to 0.90 from a depth of 8 but 1.08 to 1.31 below it past one block (six AVX2 cores, GCC 13).
constexpr std::size_t	numerical::detail::narrow_columns {8}
	The columns of C one block of gemm_rows keeps in registers over the whole depth.
constexpr std::size_t	numerical::detail::narrow_group {4}
	The rows of C one block of gemm_rows takes together: each element a chain of multiply-adds that waits on its last, four rows give the units four times the chains (M1 Pro, order 8 to 24: 1.6 to 1.9 times as fast as one row; two rows on six x86-64 AVX2 cores, the fastest there, were within 7 %).
constexpr std::size_t	numerical::detail::packed_line {64}
	The cache line a packed operand starts.
template<typename Kernel >
constexpr bool	numerical::detail::prefetches_next_row_of_c
	Whether Kernel has the C tile of the next row of tiles prefetched while a tile computes (Kernel::prefetch_rows): the tiles within a row walk one page's lines, which the processor's own prefetch follows, while the next row lands on other pages.
constexpr std::size_t	numerical::detail::split_columns {512}
	The columns of B and C, and the depth, of one real product of gemm_split: its real operands stay a few megabytes whatever the shapes.
constexpr std::size_t	numerical::detail::split_depth {256}
	The depth of one real product of gemm_split (split_columns).
constexpr std::size_t	numerical::detail::split_kept {std::size_t {1} << 18}
	The most elements an operand of gemm_split keeps from one product to the next on a thread: 2 MiB of double, every operand of a product of order 512; a larger one is let go after its product.
constexpr std::size_t	numerical::detail::split_parts_rows {64}
	The most rows of A for which gemm_split runs four real products of a block's parts on the AMX path (split_by_parts). M1 Pro, all threads, A of 16 to 64 rows by B of 256 to 1024: complex float 1.1 to 2.4 times faster, complex double 0.97 to 1.6; 96 to 128 rows by 512, complex double 1.2 to 1.4 times slower, where one wide product suited the AMX.
constexpr std::size_t	numerical::detail::plain_loop_least_width {64}
	The fewest columns of a product plain_loop_pays sends to the plain loop on the x86-64 tiles: there an order-24 SVD, whose two products of 24 × 24 × 24 it took, ran 9 % slower (GCC 13, the same bits), though such a product alone ran faster in the plain loop. On the NEON tile the small square products gain most (GCC 14, one core: a product of order 24 to 64 −32 to −50 %, a QR −5 to −21 %).
constexpr std::size_t	numerical::detail::dot_chunk {std::size_t {1} << 14U}
	The elements of one partial sum of chunked_dot: fixed, so the sum does not depend on the threads.
template<typename T >
constexpr std::size_t	numerical::detail::multiply_add_weight {1}
	The real multiply-adds one multiply-add of T weighs: 1 for a real number.
constexpr std::size_t	numerical::detail::row_group {64}
	The rows a thread takes of a matrix times a vector (matrix_times_vector).

Detailed Description

The matrix product kernel behind matrix * matrix and the einsum matmul path.

C += A·B over row-major buffers, in the loop structure of BLIS (Van Zee and van de Geijn, 2015): the depth is cut into panels of numerical::detail::gemm_kc, each operand's panel is copied once into a contiguous layout the inner loop reads at unit stride (A in slivers of numerical::detail::gemm_mr rows, B in slivers of numerical::detail::gemm_nr columns), and a register tile of gemm_mr × gemm_nr elements of C accumulates a whole panel's depth before it is written back. Each element of C is read into the tile at the start of a panel and written back at its end, so every element is summed in ascending order of the depth, exactly as the plain i-k-j loop does, on any number of threads: with the compiler's multiply-add contraction off, the two agree to the bit.

For float and double under GCC and Clang the tile is held in vectors of the target's native width (16 bytes, 32 with AVX), which both compilers keep in registers; elsewhere it is a plain array.

Definition in file gemm.hpp.

Macro Definition Documentation

NUMERICAL_OUT_OF_LINE

#define NUMERICAL_OUT_OF_LINE   [[gnu::noinline]]

Definition at line 52 of file gemm.hpp.

NUMERICAL_SCALAR_DEPTH

#define NUMERICAL_SCALAR_DEPTH

Definition at line 61 of file gemm.hpp.

NUMERICAL_UNROLL_BLOCK

#define NUMERICAL_UNROLL_BLOCK

Definition at line 69 of file gemm.hpp.