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 int)
- func ArgMin(a []float64) int
- func ArgMinAxis(dst []float64, src []float64, outer, axisLen, inner 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 Greater(dst, a, b []float64)
- func GreaterEqual(dst, a, b []float64)
- func Less(dst, a, b []float64)
- func LessEqual(dst, a, b []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 MatVecP(dst, a, v []float64, m, k int)
- func Max(a []float64) float64
- func MaxAxis(dst, src []float64, outer, axisLen, inner 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 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 int)
- 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 int)
- func SumP(a []float64) float64
- func VecMatP(dst, v, a []float64, k, n int)
- func Where(dst, cond, t, f []float64)
Constants ¶
const ( MR = 4 NR = 4 )
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.
func ArgMaxAxis ¶
ArgMaxAxis writes into dst the index (along the middle axis) of the first maximum for each [outer][inner] position. 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.
func ArgMinAxis ¶
ArgMinAxis writes into dst the index (along the middle axis) of the first minimum for each [outer][inner] position.
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 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 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 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 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].
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 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].
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 RunAxisP ¶
RunAxisP runs an axis-reduction kernel, fanning the `outer` slabs across cores above ParThreshold (measured by the total element count outer*axisLen*inner). Each worker owns a disjoint band of outer rows and writes the matching disjoint [band*inner] region of dst, reading only its own src slab — so the parallel result is identical to the serial kernel. Splitting `outer` keeps every worker's inner traversal contiguous (the cache-friendly axis), which is the case that matters: an axis-1 reduction of an (R x C) matrix has outer=R, inner=1, and the serial kernel leaves all but one core idle.
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 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 ¶
This section is empty.
Source Files