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 drop in without changing callers or tests: go-asmgen-generated kernels replace the hot ones on amd64, arm64, ppc64le and loong64 (the *_<arch> files), while these scalar versions remain the reference and the fallback elsewhere.
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 AddBlocks8F32(r *[8]float32, a []float32)
- func AddP(dst, a, b []float64)
- func ArgMax(a []float64) int
- func ArgMin(a []float64) int
- func Bin32(op int, dst, a, b []float32)
- 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 Exp(dst, src []float64)
- func ExpP(dst, src []float64)
- func HasGemm32() bool
- 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 MatMul32StridedP(dst []float32, a, b Mat32, 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 MatVecStridedP(dst, a []float64, lda int, v []float64, m, k int)
- func Max(a []float64) float64
- func Max32(a []float32) float32
- 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 Min32(a []float32) float32
- 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 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 Sqrt32(dst, src []float32)
- 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
- type Mat32
- type MatOf
Constants ¶
const ( Add32 = iota Sub32 Mul32 Div32 )
The float32 elementwise operations, in the order of the parent package's binCode: Bin32 takes these codes.
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 AddBlocks8F32 ¶ added in v0.14.0
AddBlocks8F32 adds a, a whole number of blocks of eight float32, into the eight accumulators r lane by lane: for each block i in order, r[j] += a[8i+j]. It is the inner loop of NumPy's pairwise sum (FLOAT_pairwise_sum, whose eight accumulators exist "to allow vectorization without changing summation ordering"), so a vector kernel gives exactly the same sums as the scalar loop: one IEEE float32 addition per lane per block, in the same order. It panics if len(a) is not a multiple of 8.
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 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 Bin32 ¶ added in v0.18.0
Bin32 writes a[i] OP b[i] into dst[i] for every element, OP one of Add32, Sub32, Mul32 or Div32. a and b must be at least as long as dst. A vector kernel runs where the architecture has one; each lane is one IEEE float32 operation, so the result is the scalar loop's bit for bit.
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 chunk is dot-producted with the SIMD kernel 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 HasGemm32 ¶ added in v0.26.0
func HasGemm32() bool
HasGemm32 reports whether this architecture and CPU have a float32 GEMM micro-kernel, which MatMul32StridedP requires. Without one, the caller multiplies in float64.
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 MatMul32StridedP ¶ added in v0.26.0
MatMul32StridedP is MatMulStridedP in float32, accumulating in float32 as a BLAS sgemm does. It requires HasGemm32.
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 MatVecStridedP ¶ added in v0.11.0
MatVecStridedP is MatVecP on rows lda elements apart: row i of a is a[i*lda : i*lda+k], so a sub-block of a larger row-major matrix is multiplied in place.
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 Max32 ¶ added in v0.20.0
Max32 returns the maximum of a (non-empty) as numpy.max computes it: NaN if any element is NaN, and -0 below +0 (Go's max). It reads the float32 elements directly. Converting to float64 first gives the same result (the conversion is exact) at the cost of a second array.
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 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 Sqrt32 ¶ added in v0.19.0
func Sqrt32(dst, src []float32)
Sqrt32 writes the square root of src[i] into dst[i] for every element. src must be at least as long as dst. A vector kernel runs where the architecture has one; the correctly rounded float32 square root is float32(math.Sqrt(float64(x))) bit for bit, since float64 carries more than twice float32's precision.
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 MatOf ¶ added in v0.26.0
type MatOf[T gemmFloat] struct {
Data []T
Off, RS, CS int
}
MatOf 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