Documentation
¶
Overview ¶
Package kernels holds the inner numeric loops for the ndarray package.
Every loop here is a portable, pure-Go (CGO=0) scalar kernel. They are kept behind this small, contiguous-slice API so that SIMD variants can drop in later without changing callers or tests: Phase 1 of the roadmap replaces these with go-asmgen-generated kernels across all six 64-bit Go SIMD targets (amd64, arm64, riscv64, loong64, ppc64le, s390x), selected at runtime, while these scalar versions remain the reference and the fallback.
Kernels operate on flat, contiguous []float64 slices. Shape, stride and broadcasting concerns live in the parent package; by the time a slice reaches a kernel it is already materialised in row-major order.
Index ¶
- Constants
- Variables
- func Abs(x float64) float64
- func Add(dst, a, b []float64)
- func AddP(dst, a, b []float64)
- func ArgMax(a []float64) int
- func ArgMaxAxis(dst []float64, src []float64, outer, axisLen, inner, lo, hi int)
- func ArgMin(a []float64) int
- func ArgMinAxis(dst []float64, src []float64, outer, axisLen, inner, lo, hi int)
- func Clip(dst, src []float64, lo, hi float64)
- func CumProdAxis(dst, src []float64, outer, axisLen, inner int)
- func CumSumAxis(dst, src []float64, outer, axisLen, inner int)
- func Div(dst, a, b []float64)
- func DivP(dst, a, b []float64)
- func Dot1D(a, b []float64) float64
- func Dot1DP(a, b []float64) float64
- func Equal(dst, a, b []float64)
- func Exp(dst, src []float64)
- func ExpP(dst, src []float64)
- func Greater(dst, a, b []float64)
- func GreaterEqual(dst, a, b []float64)
- func Less(dst, a, b []float64)
- func LessEqual(dst, a, b []float64)
- func Log(dst, src []float64)
- func Log10(dst, src []float64)
- func Log10P(dst, src []float64)
- func LogP(dst, src []float64)
- func Map(dst, src []float64, f func(float64) float64)
- func MapP(dst, src []float64, f func(float64) float64)
- func MatMul(dst, a, b []float64, m, k, n int)
- func MatMulP(dst, a, b []float64, m, k, n int)
- func MatMulStridedP(dst []float64, a, b Mat, m, k, n int)
- func MatVecP(dst, a, v []float64, m, k int)
- func Max(a []float64) float64
- func MaxAxis(dst, src []float64, outer, axisLen, inner, lo, hi int)
- func MaxP(a []float64) float64
- func Maximum(dst, a, b []float64)
- func Min(a []float64) float64
- func MinAxis(dst, src []float64, outer, axisLen, inner, lo, hi int)
- func MinP(a []float64) float64
- func Minimum(dst, a, b []float64)
- func Mul(dst, a, b []float64)
- func MulP(dst, a, b []float64)
- func NotEqual(dst, a, b []float64)
- func Prod(a []float64) float64
- func ProdAxis(dst, src []float64, outer, axisLen, inner, lo, hi int)
- func RepeatP(k func(dst, a, b []float64), dst, full, rep []float64, repFirst bool)
- func RunAxisP(k axisKernel, dst, src []float64, outer, axisLen, inner int)
- func SqrtP(dst, src []float64)
- func Sub(dst, a, b []float64)
- func SubP(dst, a, b []float64)
- func Sum(a []float64) float64
- func SumAxis(dst, src []float64, outer, axisLen, inner, lo, hi int)
- func SumP(a []float64) float64
- func VecMatP(dst, v, a []float64, k, n int)
- func Where(dst, cond, t, f []float64)
- type Mat
Constants ¶
const ( MR = 6 NR = 8 )
MR, NR are the micro-kernel's register-block dimensions: it computes an MR-row x NR-col tile of C per call. The packing in gemm.go matches.
const HaveReduceSIMD = true
HaveReduceSIMD reports that this build ships SIMD reduction/elementwise kernels (so the dispatch wrappers route through them rather than the scalar oracle). True on amd64/arm64; false on the four arches without vector double.
Variables ¶
var GemmThreshold = 6000
GemmThreshold is the minimum number of result elements (m*n) at which the packed GEMM fans out across goroutines. Small products run single-threaded (one worker, one pooled buffer) so the goroutine launch never dominates.
6000 is the measured crossover on the Apple-silicon bench host: below it the goroutine launch + join costs more than the parallel speedup buys, above it the fan-out wins. Concretely, 64x64 (4096 results) is ~14us serial vs ~17us parallel (serial wins) while 80x80 (6400 results) is ~26us serial vs ~22us parallel and 96x96 (9216) is ~46us serial vs ~34us parallel (parallel wins by 15-30%). The old 1<<14 (16384) floor left the whole 80x80..120x120 band — the shapes that benefit most from the cores — stuck on the slow serial path; it was sized for MatVecP (memory-bound, see matVecThreshold), not for the compute-bound GEMM. 32x32 (1024) and 64x64 (4096) stay serial, where the register-blocked single-core kernel is fastest.
var ParThreshold = 32768
ParThreshold is the minimum element count at which the parallel kernels fan out across goroutines. Below it the serial scalar kernel is used, because the goroutine launch/join cost exceeds the work. It is a var so tests can force the parallel path on small inputs (and pin it for deterministic coverage).
Functions ¶
func Abs ¶
Abs returns the absolute value of x. It exists so the parent package can route Abs through Map without importing math directly.
func AddP ¶
func AddP(dst, a, b []float64)
AddP writes a[i]+b[i] into dst[i], parallelised above ParThreshold.
func ArgMax ¶
ArgMax returns the index of the first maximum element of a (non-empty), matching numpy.argmax: ties go to the lowest index, and a NaN counts as the maximum, so the index of the first NaN wins (as Max returns NaN).
func ArgMaxAxis ¶
ArgMaxAxis writes into dst the index (along the middle axis) of the first maximum for each [outer][inner] position, or of the first NaN (see ArgMax). Layout matches the *Axis kernels.
func ArgMin ¶
ArgMin returns the index of the first minimum element of a (non-empty), matching numpy.argmin: ties go to the lowest index, and the first NaN wins.
func ArgMinAxis ¶
ArgMinAxis writes into dst the index (along the middle axis) of the first minimum for each [outer][inner] position, or of the first NaN.
func Clip ¶
Clip writes into dst each src element limited to [lo, hi], matching numpy.clip. A NaN bound or element follows Go's min/max comparison order; callers pass lo <= hi.
func CumProdAxis ¶
CumProdAxis writes the cumulative product along the middle axis into dst.
func CumSumAxis ¶
CumSumAxis writes the cumulative sum along the middle axis into dst (same shape as src), matching numpy.cumsum along an axis. Layout is [outer][axisLen][inner] as for the reductions.
func DivP ¶
func DivP(dst, a, b []float64)
DivP writes a[i]/b[i] into dst[i], parallelised above ParThreshold.
func Dot1DP ¶
Dot1DP returns the inner product of two equal-length contiguous vectors, fanned across cores above ParThreshold: each worker dot-products its chunk with the unrolled dotRange and the partials are summed. The chunked four-way grouping reassociates the sum (a few ULP, like SumP / numpy pairwise); callers needing the exact left-to-right value use the scalar Dot1D.
func Exp ¶ added in v0.2.0
func Exp(dst, src []float64)
Exp writes exp(src[i]) into dst[i]. exp1 is too large to inline, and a call per element cost half the time, so the common case (2^-54 <= |x| < 512) is exp1's main path written out in the loop; everything else calls exp1. TestExpLoopMatchesExp1 holds the two to the same bits.
func ExpP ¶ added in v0.2.0
func ExpP(dst, src []float64)
ExpP is Exp, split across cores above ParThreshold.
func Greater ¶
func Greater(dst, a, b []float64)
Greater writes the a[i]>b[i] mask (1/0) into dst[i].
func GreaterEqual ¶
func GreaterEqual(dst, a, b []float64)
GreaterEqual writes the a[i]>=b[i] mask (1/0) into dst[i].
func LessEqual ¶
func LessEqual(dst, a, b []float64)
LessEqual writes the a[i]<=b[i] mask (1/0) into dst[i].
func Log ¶ added in v0.2.0
func Log(dst, src []float64)
Log writes log(src[i]) into dst[i]. Positive normal inputs away from 1 take the table path straight from the loop; everything else goes through log1. TestLogLoopMatchesLog1 holds the two to the same bits.
func Log10 ¶ added in v0.2.0
func Log10(dst, src []float64)
Log10 writes log10(src[i]) into dst[i] as log(x)*(1/ln10), the formula of math.Log10, on this package's log.
func Log10P ¶ added in v0.2.0
func Log10P(dst, src []float64)
Log10P is Log10, split across cores above ParThreshold.
func LogP ¶ added in v0.2.0
func LogP(dst, src []float64)
LogP is Log, split across cores above ParThreshold.
func MapP ¶
MapP applies f to every element of src into dst, parallelised above ParThreshold. f must be safe for concurrent calls (the math ufuncs are).
func MatMul ¶
MatMul computes dst = a(m x k) * b(k x n) serially with the packed GEMM (the single-worker path of MatMulP). dst must be zeroed by the caller. It is kept as the named serial entry point used by callers and tests that want the exact non-parallel computation.
func MatMulP ¶
MatMulP computes dst = a(m x k) * b(k x n) with the packed, cache-blocked GEMM, fanning the m output rows across cores above GemmThreshold. dst must be zeroed by the caller.
Parallelism is the BLIS "loop-3" (ic, row-block) split with DYNAMIC, work- stealing scheduling: for every (jc, pc) cache block the B panel is packed ONCE into a shared buffer, then the m rows are diced into MR-aligned bands of MC rows and the workers pull bands off a shared atomic counter until exhausted, each packing only its own A panel and running the macro-kernel (gemmBand) over the band against the single shared packed B. Every band writes a disjoint, MR-aligned region of dst (b and the shared pb are read-only during the band), so the result is identical to the serial computation.
Two properties together gave the large-N scaling:
- Packing B once (not once per worker) removes the O(P) redundant 1-MiB repack the old row-split paid per kc block.
- Dynamic banding handles the M-series big.LITTLE asymmetry. A static one-band-per-core split makes every fast P-core wait on the slow E-cores at the join (measured: 16 cores was SLOWER than 10 on 1024², the 4 E-cores straggling). With many small bands the P-cores simply grab more of them, so adding the E-cores helps instead of hurting.
func MatMulStridedP ¶
MatMulStridedP is MatMulP over strided operand views (see Mat): dst (m x n, contiguous, zeroed by the caller) = a (m x k) * b (k x n). The views are only read, through the packing, so a transposed or sliced operand costs no copy beyond the one the packed GEMM makes anyway.
func MatVecP ¶
MatVecP computes dst = a(m x k) * v(k,) as m independent contiguous row-dots, fanned across cores above matVecThreshold. Each output dst[i] is the dot of row i of a (the contiguous slice a[i*k:(i+1)*k]) with v, so the access pattern is pure unit-stride streaming — far cheaper than routing a single-column GEMM through the packer. dst need not be pre-zeroed (each entry is written, not accumulated). Rows are disjoint, so the parallel result equals the serial one.
func Max ¶
Max returns the maximum element of a, which must be non-empty.
NaN convention: Max is NaN-propagating — if any element is NaN the result is NaN. This matches numpy.max (and numpy.maximum), Go's builtin max, and IEEE-754 maximum (not the C fmax / IEEE maxNum "ignore NaN" rule); on signed zeros it returns +0 for max(-0,+0), also matching numpy. The earlier `if v > m` form silently *ignored* NaNs (and a leading NaN poisoned the scan), diverging from numpy; this is the corrected, documented semantics.
It uses the BUILTIN max, not math.Max: both have identical NaN/signed-zero semantics, but the builtin lowers to the hardware FMAXD intrinsic on arm64 (and the equivalent elsewhere) — math.Max is an un-intrinsified call that is ~12x slower in this hot loop. The SSE2 maxSIMD kernel (MAXPD + NaN scan) is validated bit-identical to this oracle.
func MaxAxis ¶
MaxAxis reduces the middle axis by taking the maximum. Like Max it is NaN-propagating (numpy.max(axis=...)); max is exact and associative, so the inner == 1 row reduction through the SIMD max kernel is the same value.
func MaxP ¶
MaxP returns the maximum element of a (non-empty), parallelised above ParThreshold. It runs the SIMD max kernel (packed MAXPD + NaN scan on amd64; the intrinsic-lowered scalar FMAXD reducer on arm64/others) on each chunk and folds the partials with the same NaN-propagating kernel. Max is associative and NaN-propagating, so the result is exactly the serial Max — including returning NaN whenever any element is NaN (numpy.max semantics).
func Maximum ¶
func Maximum(dst, a, b []float64)
Maximum writes the pairwise maximum of a[i] and b[i] into dst[i]. It is NaN-propagating like numpy.maximum: it uses the builtin max, not math.Max, because math.Max returns +Inf for max(NaN, +Inf) where numpy returns NaN.
func Min ¶
Min returns the minimum element of a, which must be non-empty. Like Max it is NaN-propagating (any NaN -> NaN) and uses the builtin min (FMIND intrinsic), matching numpy.min including min(-0,+0) = -0.
func MinAxis ¶
MinAxis reduces the middle axis by taking the minimum, NaN-propagating (numpy.min(axis=...)).
func MinP ¶
MinP returns the minimum element of a (non-empty), parallelised above ParThreshold, via the SIMD min kernel; result identical to the serial Min, NaN-propagating like numpy.min.
func Minimum ¶
func Minimum(dst, a, b []float64)
Minimum writes the pairwise minimum of a[i] and b[i] into dst[i], NaN-propagating like numpy.minimum (builtin min; math.Min returns -Inf for min(NaN, -Inf)).
func MulP ¶
func MulP(dst, a, b []float64)
MulP writes a[i]*b[i] into dst[i], parallelised above ParThreshold.
func NotEqual ¶
func NotEqual(dst, a, b []float64)
NotEqual writes the a[i]!=b[i] mask (1/0) into dst[i].
func RepeatP ¶ added in v0.2.0
RepeatP computes dst = full OP rep, where rep is repeated along full in consecutive blocks of len(rep) — a broadcast of rep over the leading axes, such as a matrix plus a row. Each block is one call of the elementwise kernel k on live slices, so rep is never materialised to full's size; the blocks are split across cores above ParThreshold. With repFirst the operands are swapped (dst = rep OP full), for the non-commutative ops. len(dst) == len(full) is a multiple of len(rep) >= 1.
func RunAxisP ¶
RunAxisP runs an axis-reduction kernel across cores above ParThreshold (measured by the total element count outer*axisLen*inner). With at least one outer slab per worker it splits the slabs; with fewer (an axis-0 reduction has one slab) it splits the inner columns into bands of at least axisBand, which the old driver ran serially. Either way each worker writes disjoint elements of dst, each computed exactly as the serial kernel would, so the parallel result is identical to it.
func SqrtP ¶
func SqrtP(dst, src []float64)
SqrtP writes sqrt(src[i]) into dst[i], routed through the SIMD sqrt kernel (packed SQRTPD on amd64; the intrinsic-lowered scalar FSQRTD loop on arm64/others) and parallelised above ParThreshold. Unlike MapP it does NOT pass math.Sqrt through a func-pointer — that indirection blocks the compiler's FSQRTD intrinsic and the packed SSE2 kernel, the very thing that made Sqrt lose to numpy. The result is bit-identical to a scalar math.Sqrt loop.
func SubP ¶
func SubP(dst, a, b []float64)
SubP writes a[i]-b[i] into dst[i], parallelised above ParThreshold.
func SumAxis ¶
SumAxis reduces the middle axis by summation. Each row of the band is added into dst with the elementwise SIMD add, element by element in axis order, so the result is the sequential sum. With inner == 1 every output is the sum of one contiguous row, taken with the lane-parallel SIMD sum (a regrouping, as numpy's pairwise sum is): a sequential sum there is latency-bound.
func SumP ¶
SumP returns the sum of all elements of a, computed as the SIMD-summed per-chunk partials folded by a final Sum. Floating-point addition is not associative, so this grouping can differ from the strictly-sequential Sum by a few ULP — the same trade-off numpy's pairwise summation makes; callers that need the exact sequential value use Sum.
func VecMatP ¶
VecMatP computes dst = v(k,) * a(k x n) = column sums weighted by v, i.e. dst[j] = sum_p v[p]*a[p*n+j]. It accumulates row by row so the inner loop over j streams a contiguous a row and the contiguous dst, which is the cache- friendly traversal (the transpose-free analogue of mat·vec for the 1-D·2-D case). dst must be zeroed by the caller.
Types ¶
type Mat ¶
Mat is a read-only strided view of a matrix: element (i, j) is Data[Off + i*RS + j*CS]. Packing reads the operands through it, so a transposed, sliced or offset view is multiplied in place — the packing copy the GEMM makes anyway does the gather — instead of being materialised into a contiguous copy first (what BLAS's trans flags buy). Strides may be negative or zero; every (i, j) in range must index Data.
Source Files
Directories