diff --git a/Project.toml b/Project.toml index 372268b..19d280c 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" [compat] Aqua = "0.8" Bumper = "0.7" +InteractiveUtils = "1" LinearAlgebra = "1" Random = "1" SIMD = "3" @@ -26,9 +27,10 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Bumper = "8ce10254-0962-460f-a3d8-1f77fea1446e" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Random", "Strided", "Aqua", "Bumper"] +test = ["Test", "Random", "Strided", "Aqua", "Bumper", "InteractiveUtils"] diff --git a/benchmark/bench_half_packing.jl b/benchmark/bench_half_packing.jl new file mode 100644 index 0000000..7e004f8 --- /dev/null +++ b/benchmark/bench_half_packing.jl @@ -0,0 +1,242 @@ +# Empirical measurement for `execute_half_packed!` (src/execution/halfpack.jl): +# B packed as usual, A read in place when the whole M extent is exactly one +# full SIMDKernel register tile (`Qm == mr(kernel)`) -- versus the fully +# packed `execute!`. +# +# julia --project=. benchmark/bench_half_packing.jl +# +# Point: profiling (see benchmark/profile_buckets.jl and the investigation +# behind this file) found packing to be a genuine 30-38% of `execute!`'s own +# time in the 8-32-element range, most of it attributable to a redundant +# copy-then-immediately-reread round trip through a freshly packed A panel +# when A already has the layout the microkernel could read directly. This +# script sweeps M fixed at exactly `mr(kernel)` (the only shape +# `execute_half_packed!` takes the in-place path for) across a range of N/K +# from small to large, to see where -- if anywhere -- skipping A's packing +# pays off and by how much, and whether that shrinks or grows with N/K. + +# This script's first sweep (below) only ever measured the simplest eligible +# case: a plain, contiguous, dense M-by-K/K-by-N matrix pair. Two follow-up +# sweeps extend that: +# +# - "non-ramp K": A's own rows are still unit-stride (eligible), but the K +# composite as a WHOLE (both operands' maps) is not a single affine ramp, +# so `_axis_of` hands the in-place A read a `PtrScatterAxis` -- a per-K-step +# indexed lookup instead of a compile-time stride multiply. This is the +# realistic case for a genuine multi-label tensor contraction (as opposed +# to a plain two-index GEMM), and it was previously only correctness- +# tested (test/execution/test_halfpack.jl), never benchmarked. +# - "alternate kernel shape": the first sweep only ever measured M pinned to +# the DEFAULT kernel's `mr` on this machine. This sweeps a second, +# non-default menu shape too. + +include(joinpath(@__DIR__, "harness.jl")) + +using QuasiStrided: execute_half_packed!, _default_kernel, kernel_shapes, _kernel_from_shape, + affine_ramp + +const REPS_TINY = 51 +const REPS_MAIN = 21 + +# (K, N) pairs at M = mr(kernel) exactly. Doubling ladder from tiny K/N (where +# the packing-avoidance premise should matter most, relative to a near-zero +# baseline cost) up through shapes large enough that B's own packing/compute +# dominates and any A-side saving should vanish into noise. +const HALF_PACKING_KN = [ + (1, 1), (2, 2), (4, 4), (8, 8), (16, 16), (32, 32), (63, 63), + (64, 64), (128, 128), (256, 256), (512, 512), (1024, 1024), + (2048, 2048), (4096, 512), (512, 4096), +] + +reps_for(K::Int, N::Int) = (K * N <= 64 * 64) ? REPS_TINY : REPS_MAIN + +function time_packed(::Type{T}, fx; reps::Int) where {T} + plan = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC) + execute!(plan, one(T), zero(T)) + t = median_time_s(() -> execute!(plan, one(T), zero(T)); reps = reps) + return t +end + +function time_half_packed(::Type{T}, fx; reps::Int) where {T} + plan = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC) + execute_half_packed!(plan, one(T), zero(T)) + t = median_time_s(() -> execute_half_packed!(plan, one(T), zero(T)); reps = reps) + return (t, plan) +end + +# --------------------------------------------------------------------------- +# "Non-ramp K" fixture: A[m,k1,k2] dense (m fastest, then k1, then k2) so A's +# OWN map is a clean ramp; B[k2,k1,n] dense (k2 fastest) so B's map, over the +# SAME (k1,k2) label order, is not -- stride(k2 in B) = 1 is far smaller than +# stride(k1 in B) = k2n, the opposite of what folding k1,k2 into one ramp +# needs. `affine_ramp` requires EVERY map to fold, so the whole K composite +# (shared by both operands) classifies as non-ramp, even though A's own +# addresses alone would trivially ramp -- exactly the generic multi-label +# case, and the in-place A read pays for it (indexed `PtrScatterAxis` lookups, +# not a compile-time stride multiply) regardless of which operand "caused" it. +function build_nonramp_k(::Type{T}, M::Int, k1n::Int, k2n::Int, N::Int, rng) where {T} + Amat = randn(rng, T, M, k1n, k2n) + Bmat = randn(rng, T, k2n, k1n, N) + Cmat = zeros(T, M, N) + Av = StridedView(Amat) + Bv = StridedView(Bmat) + Cv = StridedView(Cmat) + # indA labels: m=1, k1=2, k2=3. indB labels B's own physical axis order + # (k2, k1, n) = (3, 2, 4). indC: (m, n) = (1, 4). + return (Av = Av, indA = (1, 2, 3), Bv = Bv, indB = (3, 2, 4), Cv = Cv, indC = (1, 4)) +end + +const NONRAMP_K1 = 4 # fixed; k2 scales to reach each target K = k1n*k2n below +const NONRAMP_TOTAL_K_N = [ + (4, 1), (8, 2), (16, 8), (32, 32), (64, 64), (128, 128), + (256, 256), (512, 512), (1024, 1024), +] + +function run_nonramp_sweep(csv, ::Type{T}, rng) where {T} + default_kernel = _default_kernel(T) + M = mr(default_kernel) + for (K, N) in NONRAMP_TOTAL_K_N + k1n = NONRAMP_K1 + k2n = cld(K, k1n) + fx = build_nonramp_k(T, M, k1n, k2n, N, rng) + plan_check = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC) + if affine_ramp(plan_check.kgroup)[1] + @warn "non-ramp K fixture unexpectedly ramps; skipping" T K N + continue + end + reps = reps_for(k1n * k2n, N) + + t_packed = time_packed(T, fx; reps = reps) + (t_half, plan) = time_half_packed(T, fx; reps = reps) + gf_packed = gflops(T, M, k1n * k2n, N, t_packed) + gf_half = gflops(T, M, k1n * k2n, N, t_half) + speedup = t_packed / t_half + + println( + csv, + "$T,$M,$(k1n * k2n),$N,$reps,$t_packed,$t_half,$gf_packed,$gf_half,$speedup,nonramp_k" + ) + @printf( + "[nonramp-K] %-10s M=%-4d K=%-5d(=%dx%d) N=%-5d packed %10.3e s (%6.2f GF/s) half %10.3e s (%6.2f GF/s) speedup %6.3fx\n", + string(T), M, k1n * k2n, k1n, k2n, N, t_packed, gf_packed, t_half, gf_half, speedup + ) + end + return nothing +end + +# --------------------------------------------------------------------------- +# Alternate (non-default) kernel shapes: the first and non-ramp sweeps above +# only ever pin M to mr() of the DEFAULT kernel this machine's ISA resolves +# to. Sweep every OTHER menu shape too (kernel_shapes(T), real method), at a +# fixed mid-range (K, N), to check the win isn't an artifact of one specific +# register shape. +const ALT_SHAPE_K, ALT_SHAPE_N = 128, 128 + +function run_alt_shape_sweep(csv, ::Type{T}, rng) where {T} + default_kernel = _default_kernel(T) + for shape in kernel_shapes(T) + kernel = _kernel_from_shape(shape, T) + M = mr(kernel) + M == mr(default_kernel) && shape == (mr(default_kernel), nr(default_kernel), lanewidth(default_kernel)) && continue + spec = ShapeSpec("alt_$(M)x$(ALT_SHAPE_K)x$(ALT_SHAPE_N)", M, ALT_SHAPE_K, ALT_SHAPE_N) + fx = build_plain(T, spec, rng) + reps = REPS_MAIN + + plan_p = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC; kernel = kernel) + execute!(plan_p, one(T), zero(T)) + t_packed = median_time_s(() -> execute!(plan_p, one(T), zero(T)); reps = reps) + plan_h = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC; kernel = kernel) + execute_half_packed!(plan_h, one(T), zero(T)) + t_half = median_time_s(() -> execute_half_packed!(plan_h, one(T), zero(T)); reps = reps) + + gf_packed = gflops(T, M, ALT_SHAPE_K, ALT_SHAPE_N, t_packed) + gf_half = gflops(T, M, ALT_SHAPE_K, ALT_SHAPE_N, t_half) + speedup = t_packed / t_half + + println( + csv, + "$T,$M,$ALT_SHAPE_K,$ALT_SHAPE_N,$reps,$t_packed,$t_half,$gf_packed,$gf_half,$speedup,alt_shape_$(shape[1])x$(shape[2])W$(shape[3])" + ) + @printf( + "[alt-shape] %-10s kernel=%dx%d/W%d packed %10.3e s (%6.2f GF/s) half %10.3e s (%6.2f GF/s) speedup %6.3fx\n", + string(T), shape[1], shape[2], shape[3], t_packed, gf_packed, t_half, gf_half, speedup + ) + end + return nothing +end + +const OUTDIR = results_dir() +mkpath(OUTDIR) +const CSV_PATH = joinpath(OUTDIR, "half_packing.csv") +const PROV_PATH = joinpath(OUTDIR, "half_packing_PROVENANCE.txt") + +function main() + print_env_header(stdout, "bench_half_packing.jl") + csv = open(CSV_PATH, "w") + println(csv, "dtype,M,K,N,reps,t_packed,t_half_packed,gf_packed,gf_half_packed,speedup_half_over_packed,fixture") + + canaries = Float64[] + rng = Random.MersenneTwister(0x4A1FBACC) + push!(canaries, run_canary(rng, "start")) + + for T in DTYPES # Float64, Float32 -- execute_half_packed! is real-only by design + # M is fixed at exactly mr(kernel) for the DEFAULT kernel this dtype + # resolves to on this machine -- the only shape the in-place path takes. + default_kernel = _default_kernel(T) + M = mr(default_kernel) + for (K, N) in HALF_PACKING_KN + spec = ShapeSpec("$(M)x$(K)x$(N)", M, K, N) + fx = build_plain(T, spec, rng) + reps = reps_for(K, N) + + t_packed = time_packed(T, fx; reps = reps) + (t_half, plan) = time_half_packed(T, fx; reps = reps) + gf_packed = gflops(T, M, K, N, t_packed) + gf_half = gflops(T, M, K, N, t_half) + speedup = t_packed / t_half + + println( + csv, + "$T,$M,$K,$N,$reps,$t_packed,$t_half,$gf_packed,$gf_half,$speedup,plain" + ) + @printf( + "%-10s M=%-4d (mr) K=%-5d N=%-5d packed %10.3e s (%6.2f GF/s) half %10.3e s (%6.2f GF/s) half/packed speedup %6.3fx kernel=%dx%d/W%d\n", + string(T), M, K, N, t_packed, gf_packed, t_half, gf_half, speedup, + mr(plan.kernel), nr(plan.kernel), lanewidth(plan.kernel) + ) + end + + run_nonramp_sweep(csv, T, rng) + run_alt_shape_sweep(csv, T, rng) + end + + push!(canaries, run_canary(rng, "end")) + close(csv) + + spread = relative_spread(canaries) + @printf("\ncanary spread: %.1f%% %s\n", 100spread, canaries) + if spread > 0.1 + @warn "canary spread exceeds 10%: the machine was not quiet enough for " * + "a few-percent conclusion. Re-run before believing any of these ratios." + end + + open(PROV_PATH, "w") do io + print_env_header(io, "bench_half_packing.jl") + println(io, "git commit = ", git_commit()) + println(io, "canaries = ", canaries) + @printf(io, "canary spread = %.2f%%\n", 100spread) + println( + io, + "\nM is fixed per dtype at mr(default_kernel(T)) on THIS machine -- see the", + ) + println( + io, + "printed kernel= shape per row. speedup_half_over_packed > 1 means the", + ) + println(io, "in-place-A path won; read half_packing.csv for the full sweep.") + end + println("\nwrote ", CSV_PATH) + return nothing +end + +main() diff --git a/benchmark/bench_skip_packing.jl b/benchmark/bench_skip_packing.jl new file mode 100644 index 0000000..828d6cd --- /dev/null +++ b/benchmark/bench_skip_packing.jl @@ -0,0 +1,133 @@ +# Empirical crossover sweep for `execute_direct!` (src/execution/direct.jl): +# the unpacked, scalar, BLIS-SUP-style path for small contractions, versus the +# existing packed `execute!`. +# +# julia --project=. benchmark/bench_skip_packing.jl +# +# Point: `plan_contract`/`execute!` always pack A/B into register-format +# panels before running the SIMD microkernel, even for tiny shapes where +# packing's O(MK+NK) cost isn't repaid by the microkernel's cache reuse. +# `execute_direct!` skips packing and the microkernel entirely (a plain scalar +# triple loop over the same AxisGroup-derived offsets `execute!` uses). This +# script sweeps shapes from tiny (M=N=K=1) up past where the packed path +# clearly wins, to find the empirical crossover -- per the project's +# "no claim without a measurement" rule, no threshold is hardcoded anywhere in +# src/ off the back of this file; a human reads the CSV/ranking and decides. +# +# Shapes deliberately include genuinely small M *and* N *and* K together (not +# just one small free extent, as `SMALL_SHAPES` in harness.jl does), plus a +# small-K family around the `dim63_2_1_2` regression shape (K=63) at a few +# M/N sizes, since that is the shape that motivated this work. + +include(joinpath(@__DIR__, "harness.jl")) + +using QuasiStrided: execute_direct! + +const REPS_TINY = 51 # small shapes need more reps: timer-resolution noise +const REPS_MAIN = 21 # matches bench_complex_efficiency.jl's REPS + +# (M, K, N) triples. Doubling ladder from 1 up to past the packed crossover, +# plus asymmetric small-K cases (the dim63_2_1_2 family) and small-M/N cases +# (mirroring SMALL_SHAPES's motivation) at a few different absolute scales. +const SKIP_PACKING_SHAPES = [ + # --- every dim tiny: doubling ladder --- + (1, 1, 1), (2, 2, 2), (4, 4, 4), (8, 8, 8), (16, 16, 16), + (32, 32, 32), (64, 64, 64), (128, 128, 128), (256, 256, 256), + # --- small-K family (dim63_2_1_2's K=63 is the motivating shape) --- + (32, 63, 32), (64, 63, 64), (128, 63, 128), (256, 63, 256), + (32, 1, 32), (64, 1, 64), (256, 1, 256), + (32, 8, 32), (64, 8, 64), (256, 8, 256), + # --- small-M or small-N, K not tiny --- + (4, 64, 256), (256, 64, 4), (4, 256, 256), (256, 256, 4), + (12, 64, 256), (256, 64, 12), + # --- small M and N together, K larger --- + (8, 256, 8), (16, 128, 16), (4, 512, 4), +] + +reps_for(M::Int, K::Int, N::Int) = (M * K * N <= 64^3) ? REPS_TINY : REPS_MAIN + +function time_packed(::Type{T}, fx; reps::Int) where {T} + plan = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC) + execute!(plan, one(T), zero(T)) + t = median_time_s(() -> execute!(plan, one(T), zero(T)); reps = reps) + return t +end + +function time_direct(::Type{T}, fx; reps::Int) where {T} + plan = plan_contract(fx.Cv, fx.Av, fx.indA, fx.Bv, fx.indB, fx.indC) + execute_direct!(plan, one(T), zero(T)) + t = median_time_s(() -> execute_direct!(plan, one(T), zero(T)); reps = reps) + return t +end + +const OUTDIR = results_dir() +mkpath(OUTDIR) +const CSV_PATH = joinpath(OUTDIR, "skip_packing.csv") +const PROV_PATH = joinpath(OUTDIR, "skip_packing_PROVENANCE.txt") + +function main() + print_env_header(stdout, "bench_skip_packing.jl") + csv = open(CSV_PATH, "w") + println(csv, "dtype,M,K,N,reps,t_packed,t_direct,gf_packed,gf_direct,speedup_direct_over_packed") + + canaries = Float64[] + rng = Random.MersenneTwister(0xDEC0DE01) + push!(canaries, run_canary(rng, "start")) + + for T in DTYPES # Float64, Float32 only -- execute_direct! is exercised + # for complex types in test/execution/test_direct.jl; this sweep is + # about the packed-vs-direct crossover, which the real kernels' + # measured menu already anchors, so keep the sweep itself real-only. + for (M, K, N) in SKIP_PACKING_SHAPES + spec = ShapeSpec("$(M)x$(K)x$(N)", M, K, N) + fx = build_plain(T, spec, rng) + reps = reps_for(M, K, N) + + t_packed = time_packed(T, fx; reps = reps) + t_direct = time_direct(T, fx; reps = reps) + gf_packed = gflops(T, M, K, N, t_packed) + gf_direct = gflops(T, M, K, N, t_direct) + speedup = t_packed / t_direct + + println( + csv, + "$T,$M,$K,$N,$reps,$t_packed,$t_direct,$gf_packed,$gf_direct,$speedup" + ) + @printf( + "%-10s %4dx%4dx%4d packed %10.3e s (%6.2f GF/s) direct %10.3e s (%6.2f GF/s) direct/packed speedup %6.3fx\n", + string(T), M, K, N, t_packed, gf_packed, t_direct, gf_direct, speedup + ) + end + end + + push!(canaries, run_canary(rng, "end")) + close(csv) + + spread = relative_spread(canaries) + @printf("\ncanary spread: %.1f%% %s\n", 100spread, canaries) + if spread > 0.1 + @warn "canary spread exceeds 10%: the machine was not quiet enough for " * + "a few-percent conclusion. Re-run before believing any crossover read " * + "off this CSV." + end + + open(PROV_PATH, "w") do io + print_env_header(io, "bench_skip_packing.jl") + println(io, "git commit = ", git_commit()) + println(io, "canaries = ", canaries) + @printf(io, "canary spread = %.2f%%\n", 100spread) + println( + io, + "\nNo threshold is derived automatically here -- read speedup_direct_over_packed", + ) + println( + io, + "in skip_packing.csv against M,K,N and pick the crossover by hand, per shape", + ) + println(io, "family (all-tiny vs small-K vs small-M/N).") + end + println("\nwrote ", CSV_PATH) + return nothing +end + +main() diff --git a/benchmark/probes/fmaddsub_codegen.jl b/benchmark/probes/fmaddsub_codegen.jl index 7b34f0d..f559d95 100644 --- a/benchmark/probes/fmaddsub_codegen.jl +++ b/benchmark/probes/fmaddsub_codegen.jl @@ -75,31 +75,41 @@ kernels(::Type{ComplexF64}) = [ ("planar", PlanarKernel(Val(24), Val(3), ComplexF64, Val(8))), ("1m", OneMKernel(Val(4), Val(6), ComplexF64, Val(4))), ("1m", OneMKernel(Val(12), Val(8), ComplexF64, Val(8))), - [("fmaddsub", FMAddSubKernel(Val(MR), Val(NR), ComplexF64, Val(W))) - for (MR, NR, W) in QuasiStrided.KERNEL_SHAPES_C64_FMADDSUB]..., + [ + ("fmaddsub", FMAddSubKernel(Val(MR), Val(NR), ComplexF64, Val(W))) + for (MR, NR, W) in QuasiStrided.KERNEL_SHAPES_C64_FMADDSUB + ]..., ] kernels(::Type{ComplexF32}) = [ ("planar", PlanarKernel(Val(8), Val(5), ComplexF32, Val(8))), ("planar", PlanarKernel(Val(48), Val(3), ComplexF32, Val(16))), ("1m", OneMKernel(Val(8), Val(6), ComplexF32, Val(8))), ("1m", OneMKernel(Val(24), Val(8), ComplexF32, Val(16))), - [("fmaddsub", FMAddSubKernel(Val(MR), Val(NR), ComplexF32, Val(W))) - for (MR, NR, W) in QuasiStrided.KERNEL_SHAPES_C32_FMADDSUB]..., + [ + ("fmaddsub", FMAddSubKernel(Val(MR), Val(NR), ComplexF32, Val(W))) + for (MR, NR, W) in QuasiStrided.KERNEL_SHAPES_C32_FMADDSUB + ]..., ] -println("cpu = ", Sys.CPU_NAME, " julia = ", VERSION, +println( + "cpu = ", Sys.CPU_NAME, " julia = ", VERSION, " (LLVM target: ", Base.JLOptions().cpu_target == C_NULL ? "native" : - unsafe_string(Base.JLOptions().cpu_target), ")") -@printf("%-10s %-9s %-10s %8s %8s %5s %5s %4s %7s %5s %5s %5s %8s %8s\n", + unsafe_string(Base.JLOptions().cpu_target), ")" +) +@printf( + "%-10s %-9s %-10s %8s %8s %5s %5s %4s %7s %5s %5s %5s %8s %8s\n", "T", "method", "shape", "fmaddsub", "fmsubadd", "fma", "fms", "mul", "add/sub", - "shuf", "blend", "bcast", "spill_st", "spill_ld") + "shuf", "blend", "bcast", "spill_st", "spill_ld" +) for T in (ComplexF64, ComplexF32) for (name, k) in kernels(T) r = probe(k) shape = "$(QuasiStrided.mr(k))x$(QuasiStrided.nr(k))/W$(QuasiStrided.lanewidth(k))" - @printf("%-10s %-9s %-10s %8d %8d %5d %5d %4d %7d %5d %5d %5d %8d %8d\n", + @printf( + "%-10s %-9s %-10s %8d %8d %5d %5d %4d %7d %5d %5d %5d %8d %8d\n", T, name, shape, r.fmaddsub, r.fmsubadd, r.fma, r.fms, r.mul, r.addsub, - r.shuf, r.blend, r.bcast, r.spill_st, r.spill_ld) + r.shuf, r.blend, r.bcast, r.spill_st, r.spill_ld + ) if DUMP println(r.asm) end diff --git a/benchmark/submit_half_packing.sh b/benchmark/submit_half_packing.sh new file mode 100755 index 0000000..c9eebd6 --- /dev/null +++ b/benchmark/submit_half_packing.sh @@ -0,0 +1,52 @@ +#!/bin/bash +# Slurm job for benchmark/bench_half_packing.jl: the packed-vs-half-packed +# (in-place A) sweep at M = mr(kernel) exactly. Single core, single node -- +# no internal parallelism, so plain sbatch is the right tool (not disBatch). +# Pick the microarchitecture with --constraint, one job per node class, since +# mr(kernel) itself is ISA-dependent (AVX2 vs AVX-512 pick different default +# register shapes) and any win/loss ratio may be too: +# +# sbatch --reservation=rocky8 --constraint=rome benchmark/submit_half_packing.sh # Zen2, AVX2 +# sbatch --reservation=rocky8 --constraint=icelake benchmark/submit_half_packing.sh # Ice Lake-SP, AVX-512 +# sbatch --reservation=rocky8 --constraint=genoa benchmark/submit_half_packing.sh # Zen4, AVX-512 +# +# Submit from the repo root. See https://wiki.flatironinstitute.org/SCC/Software/Slurm. +# Check the job's `cpu =` header and the printed kernel= shape per row before +# trusting a conclusion. +#SBATCH --job-name=qs-half-packing +#SBATCH --partition=ccq +#SBATCH --nodes=1 +#SBATCH --ntasks=1 +#SBATCH --cpus-per-task=1 +#SBATCH --mem=16G +#SBATCH --time=01:00:00 +#SBATCH --output=benchmark/results/slurm-%j.out + +set -euo pipefail +cd "${SLURM_SUBMIT_DIR:?SLURM_SUBMIT_DIR not set -- run this script via sbatch, not directly}" + +# See submit_blocking_model.sh / submit_complex_efficiency.sh / submit_skip_packing.sh +# for why this resolution block: on Rocky 8 nodes (--reservation=rocky8) the +# module system's own binaries need a newer glibc than the node has, so +# `module load` fails and there is no `julia`; fall back to $QS_JULIA or a +# juliaup-installed 1.12 on the shared home (official Julia builds need only +# glibc 2.17). +module load modules/2.5-beta1 || true +module load julia/1.12.6 || true +MODJULIA=$(command -v julia 2>/dev/null || true) +case "$MODJULIA" in /mnt/sw/*) ;; *) MODJULIA="" ;; esac # only a module-provided one +JULIA=${QS_JULIA:-${MODJULIA:-$(ls -d "$HOME"/.julia/juliaup/julia-1.12*/bin/julia 2>/dev/null | tail -1)}} +[ -x "${JULIA:-}" ] || { echo "no julia found; set QS_JULIA" >&2; exit 127; } +echo "julia = $JULIA" + +export JULIA_NUM_THREADS=1 +"$JULIA" --project=. -e 'using Pkg; Pkg.instantiate()' + +"$JULIA" --project=. benchmark/bench_half_packing.jl + +echo "=== node check ===" +echo "SLURM_JOB_NODELIST = ${SLURM_JOB_NODELIST:-?}" +lscpu | grep -i "model name\|cache" || true + +echo "=== cpu check (trust the run above only if this matches the requested --constraint) ===" +grep -m1 "cpu = " "benchmark/results/"*"/half_packing_PROVENANCE.txt" || true diff --git a/benchmark/submit_skip_packing.sh b/benchmark/submit_skip_packing.sh new file mode 100755 index 0000000..8ad5cec --- /dev/null +++ b/benchmark/submit_skip_packing.sh @@ -0,0 +1,54 @@ +#!/bin/bash +# Slurm job for benchmark/bench_skip_packing.jl: the empirical packed-vs-direct +# (skip-packing) crossover sweep for small contractions. Single core, single +# node -- no internal parallelism, so plain sbatch is the right tool (not +# disBatch). Pick the microarchitecture with --constraint on the command line, +# one job per node class, since the crossover is expected to be ISA/cache-size +# dependent, not universal: +# +# sbatch --reservation=rocky8 --constraint=rome benchmark/submit_skip_packing.sh # Zen2, AVX2 +# sbatch --reservation=rocky8 --constraint=icelake benchmark/submit_skip_packing.sh # Ice Lake-SP, AVX-512 +# sbatch --reservation=rocky8 --constraint=genoa benchmark/submit_skip_packing.sh # Zen4, AVX-512 +# +# Submit from the repo root. See https://wiki.flatironinstitute.org/SCC/Software/Slurm. +# Check the job's `cpu =` header (and the "node check" section at the bottom) +# before trusting a conclusion -- confirm it actually landed on the requested +# microarchitecture. +#SBATCH --job-name=qs-skip-packing +#SBATCH --partition=ccq +#SBATCH --nodes=1 +#SBATCH --ntasks=1 +#SBATCH --cpus-per-task=1 +#SBATCH --mem=16G +#SBATCH --time=01:00:00 +#SBATCH --output=benchmark/results/slurm-%j.out + +set -euo pipefail +# Slurm runs a spooled copy of this script; SLURM_SUBMIT_DIR is the repo root. +cd "${SLURM_SUBMIT_DIR:?SLURM_SUBMIT_DIR not set -- run this script via sbatch, not directly}" + +# On Rocky 8 nodes (--reservation=rocky8) the module system's own binaries +# need a newer glibc than the node has, so `module load` fails and there is no +# `julia`; fall back to $QS_JULIA or a juliaup-installed 1.12 on the shared +# home, since the official Julia builds need only glibc 2.17. See +# submit_blocking_model.sh / submit_complex_efficiency.sh for the same +# resolution block. +module load modules/2.5-beta1 || true +module load julia/1.12.6 || true +MODJULIA=$(command -v julia 2>/dev/null || true) +case "$MODJULIA" in /mnt/sw/*) ;; *) MODJULIA="" ;; esac # only a module-provided one +JULIA=${QS_JULIA:-${MODJULIA:-$(ls -d "$HOME"/.julia/juliaup/julia-1.12*/bin/julia 2>/dev/null | tail -1)}} +[ -x "${JULIA:-}" ] || { echo "no julia found; set QS_JULIA" >&2; exit 127; } +echo "julia = $JULIA" + +export JULIA_NUM_THREADS=1 +"$JULIA" --project=. -e 'using Pkg; Pkg.instantiate()' + +"$JULIA" --project=. benchmark/bench_skip_packing.jl + +echo "=== node check ===" +echo "SLURM_JOB_NODELIST = ${SLURM_JOB_NODELIST:-?}" +lscpu | grep -i "model name\|cache" || true + +echo "=== cpu check (trust the run above only if this matches the requested --constraint) ===" +grep -m1 "cpu = " "benchmark/results/"*"/skip_packing_PROVENANCE.txt" || true diff --git a/src/QuasiStrided.jl b/src/QuasiStrided.jl index 0486c4b..5e2e633 100644 --- a/src/QuasiStrided.jl +++ b/src/QuasiStrided.jl @@ -48,10 +48,14 @@ include("planning/blocking.jl") include("execution/workspace.jl") include("planning/plan.jl") -# --- Execution: the five-loop nest, and the tile-by-tile oracle --- +# --- Execution: the five-loop nest, the tile-by-tile oracle, the opt-in +# unpacked direct path for small contractions, and the opt-in half-packed +# path (B packed, one-tile A read in place) --- include("execution/macrokernel.jl") include("execution/execute.jl") include("execution/oracle.jl") +include("execution/direct.jl") +include("execution/halfpack.jl") # --- Integrations --- include("integrations/tensoroperations.jl") diff --git a/src/execution/direct.jl b/src/execution/direct.jl new file mode 100644 index 0000000..90a9100 --- /dev/null +++ b/src/execution/direct.jl @@ -0,0 +1,138 @@ +""" + execute_direct!(plan::ContractPlan, alpha::Number, beta::Number) + +Unexported. The "skip packing" (BLIS "SUP", skinny-unpacked-style) direct path +for small contractions: computes the same result as [`execute!`](@ref) -- +`C = alpha * sum_K atransform(A) * btransform(B) + beta * C` -- as a plain +scalar triple loop that reads `A`/`B`/`C` straight from `plan.Astorage`/ +`plan.Bstorage`/`plan.Cstorage`. No packing, no register microkernel, no cache +blocking: `plan.kernel`, `plan.blocking` and `plan.workspace` are never +touched, so it works on any plan (including one built with `oracle = false`) +and has no notion of `mr`/`nr` register tiles. For shapes this small the +O(MK + NK) packing cost of the blocked path is not repaid by register-tile +reuse, which is the whole point of skipping it. + +Deliberately opt-in: nothing in [`plan_contract`](@ref) or [`contract!`](@ref) +dispatches here; call it explicitly on a plan from `plan_contract`. Its +result agrees with [`execute!`](@ref) and [`execute_tilewise!`](@ref) up to +floating-point reassociation (the K accumulation order differs). + +Semantics match `execute!`: empty output is a no-op; empty K or `alpha == 0` +applies `beta` once to every element of `C` without reading `A`/`B`; `beta` +applies exactly once per output element; `beta == 0` never reads old `C` and +`beta == 1` adds without scaling, the same BLAS-like shortcuts as the +microkernels' stores. `atransform`/`btransform` (`identity` or `conj`) are +applied to each whole loaded element, as `pack_a!`/`pack_b!` do. + +**Not allocation-free**, unlike `execute!`/`execute_tilewise!`: each call +allocates six `Vector{Int}` offset maps (lengths `Qm`, `Qm`, `Qn`, `Qn`, +`Qk`, `Qk`). Uses ordinary bounds-checked indexing throughout. Returns +`plan.Cstorage`. + +**Measured, 2026-09-25, jobs 7109278 (ccq rome/znver2, AVX2) and 7109279 (ccq +icelake-server, AVX-512); see `benchmark/bench_skip_packing.jl` and its +`benchmark/submit_skip_packing.sh`**: this path is NEVER faster than +`execute!` at any swept shape, on either ISA. Best case is near-parity at +1x1x1/2x2x2 (0.81-0.99x); by 64^3 it is 25-40x SLOWER, and it does not improve +on the small-K family (`Qk in (1, 8, 63)` at M,N up to 256) that motivated this +work either (0.02-0.18x there). The scalar loop's own throughput is pinned +near 1-1.25 GFLOP/s regardless of shape -- it never vectorizes -- while the +packed path already reaches 30-90 GFLOP/s and costs as little as ~150ns at +1x1x1. So the packing-avoidance premise this path was built to test is FALSE +as measured: at these sizes packing's O(MK+NK) cost is not what the packed +path is paying for, and skipping it while staying scalar throws away far more +(vectorization) than it saves. **Do not wire this into `plan_contract`'s +automatic dispatch or any `contract!` opt-in kwarg** -- there is no shape +where it measurably wins. It remains useful only as a third independent +correctness check alongside `execute_tilewise!`, both exercised in +`test/execution/test_direct.jl`. A real BLIS-SUP-style win would need a +still-vectorized small-shape kernel (e.g. reading contiguous/lightly-strided +operands straight into SIMD registers without packing), which is a materially +larger, unmeasured follow-up, not a variant of this function. +""" +function execute_direct!(plan::ContractPlan{T}, alpha::Number, beta::Number) where {T} + alphaT = convert(T, alpha) + betaT = convert(T, beta) + + Qm = axis_length(plan.mgroup) + Qn = axis_length(plan.ngroup) + Qk = axis_length(plan.kgroup) + + # Empty output: no-op, nothing read or written at all. + (Qm == 0 || Qn == 0) && return plan.Cstorage + + # Map order matches `_execute_nest!`: mgroup is (A, C), ngroup is (B, C), + # kgroup is (A, B). + aOffM = Vector{Int}(undef, Qm) + cOffM = Vector{Int}(undef, Qm) + fill_offsets!((aOffM, cOffM), plan.mgroup, 0, Qm) + bOffN = Vector{Int}(undef, Qn) + cOffN = Vector{Int}(undef, Qn) + fill_offsets!((bOffN, cOffN), plan.ngroup, 0, Qn) + + C = plan.Cstorage + Cbase = plan.Cbase + + # Nothing to contract: beta-only pass, A and B never read. + if Qk == 0 || iszero(alphaT) + _direct_scale_C!(C, Cbase, cOffM, cOffN, betaT) + return C + end + + aOffK = Vector{Int}(undef, Qk) + bOffK = Vector{Int}(undef, Qk) + fill_offsets!((aOffK, bOffK), plan.kgroup, 0, Qk) + + return _direct_nest!( + C, Cbase, plan.Astorage, plan.Abase, plan.Bstorage, plan.Bbase, + aOffM, cOffM, bOffN, cOffN, aOffK, bOffK, + plan.atransform, plan.btransform, alphaT, betaT + ) +end + +# Function barrier: every argument is concretely typed here, so the triple +# loop specializes on the storage types and on the transform singletons. +function _direct_nest!( + C, Cbase::Int, A, Abase::Int, B, Bbase::Int, + aOffM::Vector{Int}, cOffM::Vector{Int}, + bOffN::Vector{Int}, cOffN::Vector{Int}, + aOffK::Vector{Int}, bOffK::Vector{Int}, + atransform::FA, btransform::FB, alphaT::T, betaT::T + ) where {FA, FB, T} + Qm = length(aOffM) + Qn = length(bOffN) + Qk = length(aOffK) + for n in 1:Qn + bn = Bbase + bOffN[n] + cn = Cbase + cOffN[n] + for m in 1:Qm + am = Abase + aOffM[m] + acc = zero(T) + for k in 1:Qk + a = atransform(A[am + aOffK[k] + 1]) + b = btransform(B[bn + bOffK[k] + 1]) + acc += a * b + end + idx = cn + cOffM[m] + 1 + C[idx] = if iszero(betaT) + alphaT * acc + elseif isone(betaT) + alphaT * acc + C[idx] + else + alphaT * acc + betaT * C[idx] + end + end + end + return C +end + +# Element-by-element `C *= beta` with `scale_tile!`'s shortcuts: `beta == 1` +# touches nothing, `beta == 0` writes zeros without reading old `C`. +function _direct_scale_C!(C, Cbase::Int, cOffM::Vector{Int}, cOffN::Vector{Int}, betaT::T) where {T} + isone(betaT) && return C + for cn in cOffN, cm in cOffM + idx = Cbase + cm + cn + 1 + C[idx] = iszero(betaT) ? zero(T) : C[idx] * betaT + end + return C +end diff --git a/src/execution/halfpack.jl b/src/execution/halfpack.jl new file mode 100644 index 0000000..c399674 --- /dev/null +++ b/src/execution/halfpack.jl @@ -0,0 +1,360 @@ +# "Half-packing": pack B exactly as `execute!` does, but read A straight out +# of its own storage when it already has the one layout the SIMD microkernel +# can consume without a copy. Opt-in only, through `execute_half_packed!`. + +""" + UnpackedAView{S,C<:Axis}(storage::S, rowbase::Int, cols::C, nrows::Int) + +Stands in for a packed A micro-panel when A is read in place: `nrows` logical +rows (the kernel's `mr`), unit-stride in storage, at `axis_length(cols)` K +steps. Logical row `i`, K step `p` lives at zero-based storage address +`rowbase + i + axis_offset(cols, p)` -- the address +[`_a_step_offset`](@ref) returns for it, which `_accumulate_step` hands +straight to `panel_vload`. So unlike a `PackedPanel`, whose `panel_vload` +offset is relative to the panel pointer, this view's offset IS the full +storage address; `rowbase` already includes the operand's own base offset. +`cols` may be affine or scattered (any `Axis`); only the row direction has to +be unit-stride, because that is what makes one `Vec{W,T}` load valid. + +Borrowed, like `PackedPanel`: `cols` may be a `PtrScatterAxis` into a +workspace buffer, so the caller keeps that buffer alive. `length` is the +packed-equivalent element count `nrows * axis_length(cols)` -- what +`_execute_tile_prologue!` compares against `packed_a_length(kernel, kc)` -- +and NOT the storage length. Not bounds-checked; see +[`execute_half_packed!`](@ref) for where the check is made. +""" +struct UnpackedAView{S, C <: Axis} + storage::S + rowbase::Int + cols::C + nrows::Int +end + +Base.length(u::UnpackedAView) = u.nrows * axis_length(u.cols) + +# The K-step address of an in-place A. Does not depend on `kernel` (the +# packed formula does); the argument is there only so both methods share one +# call shape in `_accumulate_step`. +@inline _a_step_offset(u::UnpackedAView, kernel, i::Int, p::Int) = + u.rowbase + i + axis_offset(u.cols, p) + +# Dense storage only, pinned in the signature: any other storage (or an +# element-type mismatch, e.g. complex storage under a real inner kernel) is a +# MethodError here, never a wrong read. `o` is the zero-based storage address +# (tiles.jl's "address + 1 = index" convention), unchecked -- the caller has +# validated the whole A region once. +# +# DELIBERATELY NO `panel_load` METHOD: `ScalarKernel`'s K step +# (src/microkernels/scalar.jl) addresses A through `packed_a_offset` directly, +# not through `_a_step_offset`, so an `UnpackedAView` reaching it would read +# packed-layout addresses out of unpacked storage. Leaving `panel_load` +# undefined makes that a MethodError instead of a silently wrong result. +@inline function panel_vload( + ::Type{Vec{W, T}}, u::UnpackedAView{S}, o::Int + ) where {W, T, S <: DenseVector{T}} + s = u.storage + return GC.@preserve s vload(Vec{W, T}, pointer(s) + sizeof(T) * o) +end + +# Type-only half of the eligibility test; folds to a constant per plan type. +# `_copies_unchanged` is the same value gate `_pack_a_contiguous_eligible` +# uses: for a real `T` both of the driver's transforms are the identity (and +# `plan_contract` only ever stores `identity` for a real `T` anyway). +@inline _half_pack_kernel_eligible(::SIMDKernel) = true +@inline _half_pack_kernel_eligible(::Any) = false + +@inline function _half_pack_static_eligible(plan::ContractPlan{T}) where {T} + return T <: Real && _half_pack_kernel_eligible(plan.kernel) && + plan.Astorage isa DenseVector{T} && _copies_unchanged(plan.atransform, T) +end + +# The M composite's whole-block descriptors `(dM_A, dM_C)` and the eligibility +# verdict. Only called once `Qm == mr(kernel)`, so this is the single register +# sliver. For an affine-ramp composite (the common case) they come from +# `_ramp_descriptor` in closed form, `==` to what `describe_block` would +# produce (see src/execution/macrokernel.jl), and no buffer is touched; +# otherwise they are materialized into the macro M buffers (`ws.m_buf_A`/ +# `m_buf_C`, at least `mc_eff >= mr` long). An irregular `dM_C`'s +# `PtrScatterAxis` then borrows `ws.m_buf_C`, which nothing else in the +# half-packed path writes. +@inline function _half_pack_m_descriptors!(plan::ContractPlan) + ws = plan.workspace + Qm = axis_length(plan.mgroup) + (m_ramp, m_step) = affine_ramp(plan.mgroup) + (dM_A, dM_C) = if m_ramp + (_ramp_descriptor(m_step[1], 0, Qm), _ramp_descriptor(m_step[2], 0, Qm)) + else + block_descriptors!((ws.m_buf_A, ws.m_buf_C), plan.mgroup, 0, Qm) + end + ok = dM_A.regular && _unit_stride_rows(AffineAxis(dM_A.base, dM_A.stride, dM_A.count)) + return (ok, dM_A, dM_C) +end + +""" + _half_pack_a_eligible(plan::ContractPlan) -> Bool + +Whether [`execute_half_packed!`](@ref) can read A in place for `plan`. The +same shape of test as `_pack_a_contiguous_eligible` +(src/packing/pack_contiguous.jl), lifted from one sliver to the whole +contraction: real element type, a [`SIMDKernel`](@ref), dense +(`DenseVector{T}`) A storage, a copy-unchanged A transform, and the whole M +extent equal to exactly one full register tile, `axis_length(plan.mgroup) == +mr(plan.kernel)`, whose A map is a unit-stride `AffineAxis`. All but the last +two clauses fold at compile time. On a non-ramp M composite this writes the +plan's M offset buffers (scratch space every driver refills before reading). +""" +function _half_pack_a_eligible(plan::ContractPlan) + _half_pack_static_eligible(plan) || return false + axis_length(plan.mgroup) == mr(plan.kernel) || return false + return _half_pack_m_descriptors!(plan)[1] +end + +""" + execute_half_packed!(plan::ContractPlan, alpha::Number, beta::Number) + +Unexported and opt-in: nothing in [`plan_contract`](@ref) or +[`contract!`](@ref) dispatches here. Same contract and result as +[`execute!`](@ref) on every plan; on an eligible one it skips packing A. + +When [`_half_pack_a_eligible`](@ref) holds, the whole M extent is one full +`mr`-row register sliver whose rows are unit-stride in A's (dense) storage, so +each K step of A is already `mr` consecutive elements -- exactly what one +`_pack_a_contiguous!` step would copy into the packed panel. Instead of the +copy, the microkernel reads those elements in place through an +[`UnpackedAView`](@ref), which `_accumulate_step` addresses via +`_a_step_offset` (resolved on the view's type at compile time). B is packed +exactly as `execute!` packs it, and the loop nest is `execute!`'s with the +`ic` loop and the M-sliver loop each fixed at one iteration: `jc` over N in +`nc` blocks, `pc` over K in `kc` panels, then every N-sliver's micro-tile +against the one A view. K may be scattered; B and C may have any layout. + +On any ineligible plan, and for the empty-output and `Qk == 0 || alpha == 0` +short-circuits, this is exactly `execute!(plan, alpha, beta)`. The function is +therefore always safe to call: it is never wrong, only sometimes no faster. + +Per K step the arithmetic is `execute!`'s: the same `_accumulate_step`, FMA +order, K blocking and store, on the same A values read from a different +address. Results are therefore bitwise identical to `execute!`'s (and to +`execute_tilewise!`'s) in the eligible case -- observed so in +test/execution/test_halfpack.jl, which compares with `==`. + +Bounds: as in `execute!`, storage bounds are validated once per block with +[`checked_span_bounds`](@ref) -- C once per `jc` block, B and A once per +`(jc, pc)` panel -- before any read or write, and the tile calls inside go +through `unsafe_pack_b!`/`unsafe_execute_tile!`. The A check covers exactly +the addresses the in-place vector loads read (the sliver's `mr` rows against +the panel's K columns); the `UnpackedAView` loads themselves are unchecked. +Allocation-free. Returns `plan.Cstorage`. + +**Measured, 2026-09-25, jobs 7110225 (ccq rome/znver2, AVX2) and 7110226 (ccq +icelake-server, AVX-512), after rebasing onto the analytical blocking model +and `FMAddSubKernel` (main commits 454a5ec/bc2c917/c5254aa) -- superseding an +earlier measurement (jobs 7110151/7110152) taken under the old hand-tuned +blocking constants, since this path reads `plan.blocking` and could in +principle have shifted under the new model; it did not, materially: see +`benchmark/bench_half_packing.jl` and its `benchmark/submit_half_packing.sh`**, +sweeping K/N at `M = mr(kernel)` exactly (the only eligible shape): a real but +modest win, roughly 1.01-1.08x across most of the sweep, peaking around K=N +in the 16-256 range (up to ~1.05x Float64/AVX2, ~1.05x Float64/AVX-512, +~1.03-1.08x Float32/AVX-512) and settling to ~1.0-1.02x at K,N >= ~1024 as B's +packing/compute (which scales with N*K, unlike A's fixed `mr`-row packing cost +that scales with K alone) increasingly dominates the total. One regression +pocket: on AVX2 at the tiniest K/N (1-8), this path is 6-20% SLOWER than +`execute!` -- the eligibility check and loop setup cost more than the (tiny) +packing they avoid there; AVX-512 mostly does not show this (one noise-level +dip at Float64 K=N=4, ~0.93x). Net: worth keeping as opt-in exactly as +shipped; auto-wiring into `plan_contract`'s dispatch would +need an additional lower K/N gate to avoid the AVX2 regression pocket, which +has not been attempted here. + +The measurement above used the SIMPLEST eligible fixture only: a plain, +contiguous, dense M-by-K/K-by-N matrix pair at the default kernel shape. +**Extended, 2026-09-26, jobs 7111552 (rome/AVX2) and 7111553 +(icelake-server/AVX-512), same canary discipline**, to check this +generalizes, since most real multi-label tensor contractions do NOT look +like a plain two-index GEMM: + + * *Non-ramp K*: A's own rows are still unit-stride (eligible), but the K + composite as a whole (both operands' maps together) is not one affine + ramp -- e.g. `A[m,k1,k2]` dense against `B[k2,k1,n]` dense, same K label + order, incompatible per-operand strides -- so the in-place A read takes + the indexed `PtrScatterAxis` path, not a compile-time affine stride. Same + band, 1.01-1.05x typical in the mid-range, same near-parity/slight-loss + pocket at the tiniest K/N on AVX2 (0.93-0.96x); one noisy outlier + (1.46x at a tiny, sub-microsecond K=16/N=8 point on AVX-512) is a + small-sample artifact, not a new regime. + * *Alternate register shapes*: every non-default menu shape + (`kernel_shapes(T)`) also shows a positive 1.01-1.07x win at a + fixed mid-range (K,N) = (128,128) -- not an artifact of one specific + `(mr,nr,W)`. + +So the modest-win/AVX2-tiny-K-regression verdict above is not an artifact of +the simplest fixture; it holds across the realistic non-ramp-K case and +every measured register shape. +""" +function execute_half_packed!(plan::ContractPlan{T}, alpha::Number, beta::Number) where {T} + alphaT = convert(T, alpha) + betaT = convert(T, beta) + + kernel = plan.kernel + Qm = axis_length(plan.mgroup) + Qn = axis_length(plan.ngroup) + Qk = axis_length(plan.kgroup) + + # Empty output, and the beta-only pass: nothing for A to be skipped in. + (Qm == 0 || Qn == 0 || Qk == 0 || iszero(alphaT)) && return execute!(plan, alpha, beta) + + (_half_pack_static_eligible(plan) && Qm == mr(kernel)) || return execute!(plan, alpha, beta) + (ok, dM_A, dM_C) = _half_pack_m_descriptors!(plan) + ok || return execute!(plan, alpha, beta) + + ws = plan.workspace + Astorage = plan.Astorage + # `ws`: the packed B panel and every `PtrScatterAxis` borrow its buffers. + # `Astorage`: the in-place A loads are raw-pointer loads. + return GC.@preserve ws Astorage begin + _execute_half_packed_nest!( + plan, ws, kernel, dM_A, dM_C, Qn, Qk, + plan.blocking.kc, plan.blocking.nc, alphaT, betaT + ) + end +end + +# `_execute_nest!` (src/execution/execute.jl) with the M side fixed at one full +# sliver and no A packing. The B half -- descriptors, the hoisted B check, the +# pack -- is `_execute_nest!`'s verbatim. +function _execute_half_packed_nest!( + plan::ContractPlan{T}, ws, kernel::K, dM_A::BlockDescriptor, dM_C::BlockDescriptor, + Qn::Int, Qk::Int, kc_eff::Int, nc_eff::Int, alphaT::T, betaT::T + ) where {T, K} + MRk = mr(kernel) + NRk = nr(kernel) + NRp = packed_b_per_k(kernel) + + btransform = plan.btransform + + (n_ramp, n_step) = affine_ramp(plan.ngroup) + (k_ramp, k_step) = affine_ramp(plan.kgroup) + + lenA = length(plan.Astorage) + lenB = length(plan.Bstorage) + lenC = length(plan.Cstorage) + + # The one M sliver: A's rows are `dM_A`, unit-stride (checked by the + # caller); C's rows are `dM_C`, any layout. Both ranges are fixed for the + # whole call. + rng_mA = descriptor_offset_range(dM_A, ws.m_buf_A, 0) + rng_mC = descriptor_offset_range(dM_C, ws.m_buf_C, 0) + arowbase = plan.Abase + dM_A.base + + # --- loop 5: jc over N in steps of nc_eff --- + jc = 0 + while jc < Qn + nblock = min(nc_eff, Qn - jc) + n_slivers = cld(nblock, NRk) + (rng_nB, rng_nC) = if n_ramp + _ramp_slivers!( + ws.n_desc_B, ws.n_desc_C, n_step[1], n_step[2], jc, + nblock, NRk, n_slivers + ) + else + fill_offsets!((ws.n_buf_B, ws.n_buf_C), plan.ngroup, jc, nblock) + _classify_slivers!( + ws.n_desc_B, ws.n_desc_C, ws.n_buf_B, ws.n_buf_C, + nblock, NRk, n_slivers + ) + end + + # HOISTED CHECK (C) -- every micro-tile of this jc block: the one M + # sliver's rows against the union of the N slivers' columns. Its + # region does not depend on pc, so once per jc, before any write. + checked_span_bounds(plan.Cbase, rng_mC, rng_nC, lenC) + + # --- loop 4: pc over K in steps of kc_eff --- + pc = 0 + firstpanel = true + while pc < Qk + kblock = min(kc_eff, Qk - pc) + dK_A, dK_B, rng_kA, rng_kB = if k_ramp + ( + _ramp_descriptor(k_step[1], pc, kblock), + _ramp_descriptor(k_step[2], pc, kblock), + _ramp_offset_range(k_step[1], pc, kblock), + _ramp_offset_range(k_step[2], pc, kblock), + ) + else + fill_offsets!((ws.k_buf_A, ws.k_buf_B), plan.kgroup, pc, kblock) + dA = describe_block(ws.k_buf_A, 0, kblock) + dB = describe_block(ws.k_buf_B, 0, kblock) + ( + dA, dB, + descriptor_offset_range(dA, ws.k_buf_A, 0), + descriptor_offset_range(dB, ws.k_buf_B, 0), + ) + end + colsA_k = _axis_of(dK_A, ws.k_buf_A, 0) + rowsB_k = _axis_of(dK_B, ws.k_buf_B, 0) + + # HOISTED CHECK (B) -- the whole B panel of this (jc, pc), exactly + # as `_execute_nest!`'s check 1. + checked_span_bounds(plan.Bbase, rng_kB, rng_nB, lenB) + + # HOISTED CHECK (A) -- the A region this panel reads in place: the + # sliver's `mr` rows against this panel's K columns. This is the + # ONLY check covering the `UnpackedAView` vector loads, which are + # raw-pointer reads; it is exactly `_execute_nest!`'s check 2. + checked_span_bounds(plan.Abase, rng_mA, rng_kA, lenA) + + beta_eff = firstpanel ? betaT : one(T) + + # Pack the whole B panel for this (jc, pc): every N-sliver. + for s in 0:(n_slivers - 1) + sfirst = s * NRk + colsB = _axis_of(ws.n_desc_B[s + 1], ws.n_buf_B, sfirst) + bpanel = _sliver_panel(ws.packed_b, NRp, kblock, s) + _pack_sliver!( + unsafe_pack_b!, bpanel, plan.Bstorage, plan.Bbase, rowsB_k, colsB, + kernel, btransform + ) + end + + # --- loop 2: jr over N-slivers; loop 1 (ir) is the one M sliver --- + _half_packed_micro_tiles!( + kernel, plan.Cstorage, plan.Cbase, dM_C, ws, + plan.Astorage, arowbase, colsA_k, MRk, NRk, NRp, + n_slivers, kblock, alphaT, beta_eff + ) + + firstpanel = false + pc += kblock + end + + jc += nblock + end + + return plan.Cstorage +end + +# Function barrier over the K axis type (the `_axis_of` Union dies here, per +# the guardrail at the top of src/execution/macrokernel.jl), so the +# `UnpackedAView` built below is concretely typed and the micro-tile calls +# specialize on it. `unsafe_execute_micro_tile!`'s precondition is the caller's +# hoisted C check. +@inline function _half_packed_micro_tiles!( + kernel::K, Cstorage::SC, Cbase::Int, dM_C::BlockDescriptor, ws, + Astorage::SA, arowbase::Int, colsA_k::CA, MRk::Int, NRk::Int, NRp::Int, + n_slivers::Int, kblock::Int, alphaT, beta_eff + ) where {K, SC, SA, CA <: Axis} + apanel = UnpackedAView(Astorage, arowbase, colsA_k, MRk) + for s in 0:(n_slivers - 1) + sfirst = s * NRk + rowsC = _axis_of(dM_C, ws.m_buf_C, 0) + colsC = _axis_of(ws.n_desc_C[s + 1], ws.n_buf_C, sfirst) + bpanel = _sliver_panel(ws.packed_b, NRp, kblock, s) + unsafe_execute_micro_tile!( + kernel, Cstorage, Cbase, rowsC, colsC, + apanel, bpanel, kblock, alphaT, beta_eff + ) + end + return nothing +end diff --git a/src/microkernels/interface.jl b/src/microkernels/interface.jl index 0203d95..05f072a 100644 --- a/src/microkernels/interface.jl +++ b/src/microkernels/interface.jl @@ -348,9 +348,11 @@ Everything else `_execute_tile_prologue!` validates is still validated: the destination extent against the kernel's `(MR, NR)`, `kc >= 0`, both packed capacities, and the empty / `kc == 0` / `alpha == 0` short-circuits. -The one caller in this package is `_execute_nest!` (src/execution/execute.jl), which +The callers in this package are `_execute_nest!` (src/execution/execute.jl), which validates the union of an entire (ic, jc) macro block's micro-tiles in one -[`checked_span_bounds`](@ref) call before running any of them. That is an +[`checked_span_bounds`](@ref) call before running any of them, and +`_execute_half_packed_nest!` (src/execution/halfpack.jl), which does the same +per jc block for its single M sliver. That is an exactly equivalent test: the block's micro-tiles are the full cross product of its M-sliver row sets and N-sliver column sets, those sets partition the block's two offset buffers, and the check only compares range extremes -- so diff --git a/src/microkernels/simd.jl b/src/microkernels/simd.jl index a4e905d..87fdd0a 100644 --- a/src/microkernels/simd.jl +++ b/src/microkernels/simd.jl @@ -66,6 +66,18 @@ function zero_accumulator(kernel::SIMDKernel{MR, NR, T, W}) where {MR, NR, T, W} return ntuple(_ -> z, Val((MR ÷ W) * NR)) end +# Address of A's logical row `i` at K step `p`, as `_accumulate_step` hands it +# to `panel_vload`. The default -- every `packed_a` that is an actual packed +# panel (`PackedPanel`, or a plain `AbstractVector` in packed layout) -- is +# exactly the kernel's packed formula `packed_a_offset(kernel, i, p)`, so the +# packed paths (`execute!`, `execute_tilewise!`, 1m's inner kernel) compile to +# the same code as before this indirection existed. The one other method is +# `UnpackedAView`'s (src/execution/halfpack.jl), for an A read in place from +# its own storage. Dispatch is on `packed_a`'s TYPE, which `_accumulate_step` +# is specialized on anyway, so the choice is made at compile time: there is +# no runtime branch on any path. +@inline _a_step_offset(::PA, kernel, i::Int, p::Int) where {PA} = packed_a_offset(kernel, i, p) + # Fully unrolled, closure-free K-step body: one vector load per A row-vector, # one scalar load per B column, NVECA*NR FMAs, generated as straight-line code. @generated function _accumulate_step( @@ -85,7 +97,7 @@ end bvars = [Symbol(:b, j) for j in 0:(NR - 1)] load_a = [ - :($(avars[v + 1]) = panel_vload(Vec{$W, $T}, packed_a, packed_a_offset(kernel, $(v * W), p))) + :($(avars[v + 1]) = panel_vload(Vec{$W, $T}, packed_a, _a_step_offset(packed_a, kernel, $(v * W), p))) for v in 0:(NVECA - 1) ] load_b = [ @@ -117,7 +129,9 @@ is a no-op), but **not bitwise identical** (FMA grouping/order differ — compare with a tolerance, never `==`). `packed_a`/`packed_b` may be a [`PackedPanel`](@ref) — what the driver passes, and the only form that keeps a large accumulator register-resident — or any contiguous `AbstractVector{T}` -such as a `Vector` or unit-range `view`. +such as a `Vector` or unit-range `view`. `packed_a` may also be an +[`UnpackedAView`](@ref) (A read in place by `execute_half_packed!`), which +`_a_step_offset` addresses by its own formula instead of the packed one. """ @inline function Base.accumulate( kernel::SIMDKernel{MR, NR, T, W}, acc::NTuple{NV, Vec{W, T}}, diff --git a/test/execution/test_direct.jl b/test/execution/test_direct.jl new file mode 100644 index 0000000..a4ebfdb --- /dev/null +++ b/test/execution/test_direct.jl @@ -0,0 +1,343 @@ +# execute_direct!: the opt-in unpacked scalar path (src/execution/direct.jl). +# Every case is checked against BOTH execute! and execute_tilewise! (each run +# on its own freshly built plan and operands, so no path can see another's +# output or buffers) and, where the operands are plain arrays, against an +# independent label-driven reference loop. Comparisons are `≈`, never `==`: +# the three paths accumulate K in different orders. + +# Independent reference: C[indC] = alpha * sum_K opA(A[indA]) * opB(B[indB]) + +# beta * Cstart[indC], by brute force over every label assignment. Knows +# nothing about AxisGroups, strides or plans. +function _direct_ref( + A::AbstractArray, indA, B::AbstractArray, indB, Cstart::AbstractArray, indC, + alpha, beta; conjA::Bool = false, conjB::Bool = false + ) + dims = Dict{Int, Int}() + for (lbl, L) in zip(indA, size(A)) + dims[lbl] = L + end + for (lbl, L) in zip(indB, size(B)) + dims[lbl] = L + end + klabels = Tuple(l for l in indA if l in indB && !(l in indC)) + ksizes = Tuple(dims[l] for l in klabels) + T = eltype(Cstart) + out = similar(Cstart, T) + for Ic in CartesianIndices(Cstart) + at = Dict{Int, Int}(zip(indC, Tuple(Ic))) + acc = zero(T) + for Ik in CartesianIndices(ksizes) + for (l, i) in zip(klabels, Tuple(Ik)) + at[l] = i + end + a = A[(at[l] for l in indA)...] + b = B[(at[l] for l in indB)...] + acc += (conjA ? conj(a) : a) * (conjB ? conj(b) : b) + end + out[Ic] = iszero(beta) ? alpha * acc : alpha * acc + beta * Cstart[Ic] + end + return out +end + +# Run each of execute!/execute_tilewise!/execute_direct! on its own fresh +# fixture from `mk()` (which must return identical contents on every call) and +# return the three resulting C arrays. `mk()` returns +# (Cv, Av, indA, Bv, indB, indC); `plankw` goes to plan_contract verbatim. +function _direct_three_way(mk, alpha, beta; plankw...) + runners = (execute!, execute_tilewise!, execute_direct!) + return map(runners) do run! + Cv, Av, indA, Bv, indB, indC = mk() + plan = plan_contract(Cv, Av, indA, Bv, indB, indC; plankw...) + ret = run!(plan, alpha, beta) + @test ret === plan.Cstorage + Array(Cv) + end +end + +# Plain column-major matmul fixture C[m,n] = A[m,k] * B[k,n], seeded so every +# call reproduces the same data. `Cfill === nothing` gives a random start C. +function _direct_mm_maker(::Type{T}, Ma, Ka, Na, seed; Cfill = nothing) where {T} + return function () + rng = MersenneTwister(seed) + Amat = randn(rng, T, Ma, Ka) + Bmat = randn(rng, T, Ka, Na) + Cmat = Cfill === nothing ? randn(rng, T, Ma, Na) : fill(convert(T, Cfill), Ma, Na) + return (StridedView(Cmat), StridedView(Amat), (1, 2), StridedView(Bmat), (2, 3), (1, 3)) + end +end + +const _DIRECT_TYPES = (Float64, ComplexF64) +# (Ma, Ka, Na): tiny, sub-register-tile, around one tile, and several tiles +# (default Float64 kernel is at least 8x6-ish and ComplexF64's 16x6-ish on +# AVX2/AVX-512, but nothing here depends on the exact shape). +const _DIRECT_SHAPES = ( + (1, 1, 1), (1, 3, 5), (5, 3, 1), (3, 1, 2), (2, 7, 3), (7, 5, 4), + (8, 6, 6), (9, 13, 7), (17, 11, 13), (33, 40, 29), (63, 63, 2), +) + +@testset "execute_direct!: contiguous matmul shapes vs execute!/execute_tilewise!/reference ($T)" for T in _DIRECT_TYPES + ab = T <: Complex ? ((1.0, 0.0), (1.5 - 0.5im, 0.75 + 0.25im), (2.0, 1.0), (-0.5, -1.25)) : + ((1.0, 0.0), (2.5, -0.75), (1.0, 1.0), (-0.5, 1.25)) + for (idx, (Ma, Ka, Na)) in enumerate(_DIRECT_SHAPES), (alpha, beta) in ab + mk = _direct_mm_maker(T, Ma, Ka, Na, 5000 + idx) + C_exec, C_tw, C_dir = _direct_three_way(mk, alpha, beta) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + Cv0, Av0, _, Bv0, _, _ = mk() + @test C_dir ≈ _direct_ref(Array(Av0), (1, 2), Array(Bv0), (2, 3), Array(Cv0), (1, 3), alpha, beta) + end +end + +@testset "execute_direct!: M and N below one register tile of the plan's kernel ($T)" for T in _DIRECT_TYPES + # The hardcoded literals this testset used to have (e.g. Ma=5) assumed an + # AVX2/AVX-512 planar register tile size and broke on hosts where kernel + # selection picks something else at these extents -- observed on GitHub's + # hosted CI runners (a smaller fitted planar shape on some, and on + # AVX-512 hosts, main's small-M FMAddSub routing -- kernel_selection.jl's + # `_small_m_shape` -- picking a DIFFERENT kernel type at small Qm, not + # merely a smaller mr of the same one). Predicting which kernel type/shape + # kernel selection lands on for a given (Qm, Qn) is exactly the logic this + # test must not duplicate, so instead of predicting it, build the ACTUAL + # plan for each candidate shape and use it to FILTER: only candidates + # whose plan happens to land below one tile in both directions (whatever + # kernel that turns out to be, on whatever host/ISA this runs on) are + # exercised. This holds by construction on any ISA/kernel-selection + # version; the final count assertion just guards against the candidate + # list becoming vacuous (e.g. every host resolving a 1-row kernel). + candidates = ( + (1, 4, 1), (3, 5, 2), (1, 9, 5), (5, 2, 1), + (2, 4, 2), (4, 3, 1), (1, 6, 3), (6, 4, 4), + ) + exercised = 0 + for (Ma, Ka, Na) in candidates + mk = _direct_mm_maker(T, Ma, Ka, Na, 77 + Ma + 10Na) + Cv, Av, indA, Bv, indB, indC = mk() + plan = plan_contract(Cv, Av, indA, Bv, indB, indC) + (axis_length(plan.mgroup) < mr(plan.kernel) && axis_length(plan.ngroup) < nr(plan.kernel)) || continue + exercised += 1 + C_exec, C_tw, C_dir = _direct_three_way(mk, 1.25, 0.5) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + end + @test exercised > 0 +end + +@testset "execute_direct!: ignores kernel shape and works without oracle buffers" begin + Random.seed!(424242) + Ma, Ka, Na = 11, 13, 7 + Amat, Bmat, Cstart = randn(Ma, Ka), randn(Ka, Na), randn(Ma, Na) + expected = 2.0 .* (Amat * Bmat) .+ 0.5 .* Cstart + for kw in ( + (kernel = ScalarKernel(Val(4), Val(3), Float64), kc = 5), + (kernel = SIMDKernel(Val(8), Val(6), Float64),), + (oracle = false,), + ) + Cmat = copy(Cstart) + plan = _mm_plan(Cmat, Amat, Bmat; kw...) + # The direct path must not touch the packing workspace at all. + # `isequal`, not `==`: the buffers are `undef`, so their (reused-heap) + # contents may hold NaN bit patterns, and `NaN == NaN` is false. + pa, pb = copy(plan.workspace.packed_a), copy(plan.workspace.packed_b) + execute_direct!(plan, 2.0, 0.5) + @test Cmat ≈ expected + @test isequal(plan.workspace.packed_a, pa) + @test isequal(plan.workspace.packed_b, pb) + end +end + +@testset "execute_direct!: K = 0 applies beta once, never reads A/B ($T)" for T in _DIRECT_TYPES + Ma, Na = 5, 4 + for beta in (0.5, 1.0, 0.0, -2.0) + Cstart = randn(MersenneTwister(31), T, Ma, Na) + mk = () -> ( + StridedView(copy(Cstart)), StridedView(fill(convert(T, NaN), Ma, 0)), (1, 2), + StridedView(fill(convert(T, Inf), 0, Na)), (2, 3), (1, 3), + ) + C_exec, C_tw, C_dir = _direct_three_way(mk, 1.0, beta) + @test C_dir ≈ beta .* Cstart + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + @test all(isfinite, C_dir) + end +end + +@testset "execute_direct!: alpha = 0 applies beta once, never reads A/B ($T)" for T in _DIRECT_TYPES + Ma, Ka, Na = 5, 6, 4 + for beta in (1.25, 1.0, 0.0) + Cstart = randn(MersenneTwister(32), T, Ma, Na) + mk = () -> ( + StridedView(copy(Cstart)), StridedView(fill(convert(T, NaN), Ma, Ka)), (1, 2), + StridedView(fill(convert(T, Inf), Ka, Na)), (2, 3), (1, 3), + ) + C_exec, C_tw, C_dir = _direct_three_way(mk, 0.0, beta) + @test C_dir ≈ beta .* Cstart + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + @test all(isfinite, C_dir) + end +end + +@testset "execute_direct!: beta = 0 never reads old C, like execute! ($T)" for T in _DIRECT_TYPES + # The engine's convention (src/microkernels/interface.jl `_axpby_tile!`, + # `scale_tile!`): `beta == 0` overwrites without loading, so a NaN-poisoned + # C must not leak into the result on any of the three paths. + for (Ma, Ka, Na) in ((1, 1, 1), (7, 5, 4), (17, 11, 13)) + mk = _direct_mm_maker(T, Ma, Ka, Na, 900 + Ma; Cfill = NaN) + C_exec, C_tw, C_dir = _direct_three_way(mk, 1.5, 0.0) + @test all(isfinite, C_dir) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + end + # ... and the alpha = 0 / beta = 0 beta-only pass writes exact zeros. + Cv, Av, indA, Bv, indB, indC = _direct_mm_maker(T, 4, 3, 5, 901; Cfill = NaN)() + execute_direct!(plan_contract(Cv, Av, indA, Bv, indB, indC), 0.0, 0.0) + @test all(iszero, Array(Cv)) +end + +@testset "execute_direct!: empty output is a no-op ($T)" for T in _DIRECT_TYPES + Cmat = fill(convert(T, NaN), 3, 0) + plan = plan_contract( + StridedView(Cmat), StridedView(randn(T, 3, 4)), (1, 2), + StridedView(randn(T, 4, 0)), (2, 3), (1, 3) + ) + @test execute_direct!(plan, 1.0, 2.0) === plan.Cstorage + @test size(Cmat) == (3, 0) + + Cmat2 = fill(convert(T, NaN), 0, 3) + plan2 = plan_contract( + StridedView(Cmat2), StridedView(randn(T, 0, 4)), (1, 2), + StridedView(randn(T, 4, 3)), (2, 3), (1, 3) + ) + @test execute_direct!(plan2, 1.0, 2.0) === plan2.Cstorage + @test size(Cmat2) == (0, 3) +end + +@testset "execute_direct!: worked fixture, plain and permuted views" begin + A, B, Cref = _worked_fixture() + C = zeros(Float64, 3, 4, 2) + plan = plan_contract(StridedView(C), StridedView(A), _INDA, StridedView(B), _INDB, _INDC) + execute_direct!(plan, 1.0, 0.0) + @test C ≈ Cref + + # Mirrors test_execute.jl's permuted-view case: Ap[k,b,a] == A[a,k,b], + # Bp[n,k] == B[k,n]. + Ap = permutedims(A, (2, 3, 1)) + Bp = permutedims(B, (2, 1)) + C2 = zeros(Float64, 3, 4, 2) + plan2 = plan_contract(StridedView(C2), StridedView(Ap), (2, 3, 1), StridedView(Bp), (4, 2), _INDC) + execute_direct!(plan2, 1.0, 0.0) + @test C2 ≈ Cref + + # Sliced (nonzero-offset) A and C, mirroring the sliced-destination case. + A3 = reshape(collect(1.0:40.0), 4, 5, 2) + Cfull = zeros(Float64, 4, 4, 2) + Av = StridedView(view(A3, 2:4, 1:5, 1:2)) + Cv = StridedView(view(Cfull, 2:4, :, :)) + @test offset(Av) != 0 && offset(Cv) != 0 + plan3 = plan_contract(Cv, Av, _INDA, StridedView(B), _INDB, _INDC) + execute_direct!(plan3, 1.0, 0.0) + Cref3 = [sum(A3[a, k, b] * B[k, n] for k in 1:5) for a in 2:4, n in 1:4, b in 1:2] + @test Array(Cv) ≈ Cref3 + @test all(iszero, Cfull[1, :, :]) # row outside the slice untouched +end + +@testset "execute_direct!: permuted A, negative-stride B, sliced-offset C ($T)" for T in _DIRECT_TYPES + # helpers.jl's `scattered_fixture` (benchmark/harness.jl's `build_scattered` + # shape): A[a,k,b] with a zero stride on b, permuted to (k,b,a); B[k,n] + # with a negative k stride; C[a,n,b] a strided slice of a padded array. + sentinel = convert(T, 7) + for (a_n, k_n, b_n, n_n) in ((32, 32, 8, 32), (3, 5, 2, 4), (1, 1, 1, 1), (9, 2, 3, 1)) + for (alpha, beta) in ((1.0, 0.0), (-1.5, 0.5)) + Cstart = randn(MersenneTwister(40 + a_n), T, a_n, n_n, b_n) + mk = function () + Cv, Av, indA, Bv, indB, indC = scattered_fixture(T, a_n, k_n, b_n, n_n) + fill!(parent(Cv), sentinel) + copyto!(Cv, Cstart) + return (Cv, Av, indA, Bv, indB, indC) + end + C_exec, C_tw, C_dir = _direct_three_way(mk, alpha, beta) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + Cv, Av, indA, Bv, indB, indC = mk() + @test C_dir ≈ _direct_ref(Array(Av), indA, Array(Bv), indB, Cstart, indC, alpha, beta) + + # Nothing outside C's slice is written: reset the slice to the + # sentinel after a direct run and the whole parent must be uniform. + plan = plan_contract(Cv, Av, indA, Bv, indB, indC) + @test offset(Cv) != 0 + # (StridedViews canonicalizes a length-1 axis's stride, so the + # negative stride only survives when k_n > 1.) + @test any(<(0), strides(Bv)) || k_n == 1 + execute_direct!(plan, alpha, beta) + @test Array(Cv) ≈ C_dir + fill!(Cv, sentinel) + @test all(==(sentinel), parent(Cv)) + end + end +end + +@testset "execute_direct!: operand-role swap (C stored N-major) ($T)" for T in _DIRECT_TYPES + # C[n,m] column-major: the M composite is not unit-stride in C but the N + # one is, so for a real eltype plan_contract swaps roles and + # plan.Astorage is B's parent. The swap is real-only (plan.jl), so the + # complex case instead covers an M composite that is strided in C. + Ma, Ka, Na = 5, 7, 19 + mk = function () + rng = MersenneTwister(606) + Amat = randn(rng, T, Ma, Ka) + Bmat = randn(rng, T, Ka, Na) + Cmat = randn(rng, T, Na, Ma) + return (StridedView(Cmat), StridedView(Amat), (1, 2), StridedView(Bmat), (2, 3), (3, 1)) + end + Cv, Av, indA, Bv, indB, indC = mk() + plan = plan_contract(Cv, Av, indA, Bv, indB, indC) + # Premise: the swap happened exactly when the planner says it should. + @test plan.Astorage === (T <: Real ? parent(Bv) : parent(Av)) + for (alpha, beta) in ((1.0, 0.0), (0.5, -1.0)) + C_exec, C_tw, C_dir = _direct_three_way(mk, alpha, beta) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + @test C_dir ≈ _direct_ref(Array(Av), indA, Array(Bv), indB, Array(Cv), indC, alpha, beta) + end +end + +@testset "execute_direct!: conjA/conjB, flags and conj views ($T)" for T in (ComplexF64, ComplexF32) + for (Ma, Ka, Na) in ((1, 1, 1), (3, 4, 2), (17, 9, 13)), + conjA in (false, true), conjB in (false, true) + alpha, beta = (0.75 - 1.0im), (0.5 + 0.5im) + mk = _direct_mm_maker(T, Ma, Ka, Na, 3000 + Ma) + C_exec, C_tw, C_dir = _direct_three_way(mk, alpha, beta; conjA = conjA, conjB = conjB) + @test C_dir ≈ C_exec + @test C_dir ≈ C_tw + Cv, Av, _, Bv, _, _ = mk() + @test C_dir ≈ _direct_ref( + Array(Av), (1, 2), Array(Bv), (2, 3), Array(Cv), (1, 3), alpha, beta; + conjA = conjA, conjB = conjB + ) + # A conj-wrapped view composes with the flag by XOR; the plan's + # transform must come out the same, and so must the result. + mkv = function () + Cv, Av, iA, Bv, iB, iC = mk() + return (Cv, conj(Av), iA, Bv, iB, iC) + end + Cvv, Avv, iA, Bvv, iB, iC = mkv() + planv = plan_contract(Cvv, Avv, iA, Bvv, iB, iC; conjA = !conjA, conjB = conjB) + execute_direct!(planv, alpha, beta) + @test Array(Cvv) ≈ C_dir + end +end + +@testset "execute_direct!: repeated calls on one plan, real-typed alpha into complex" begin + # A plan is reusable: a second call accumulates onto the first's output. + T = ComplexF64 + mk = _direct_mm_maker(T, 6, 5, 4, 1234) + Cv, Av, indA, Bv, indB, indC = mk() + Cstart = Array(Cv) + plan = plan_contract(Cv, Av, indA, Bv, indB, indC) + execute_direct!(plan, 2, 1) + execute_direct!(plan, 2, 1) + AB = Array(Av) * Array(Bv) + @test Array(Cv) ≈ Cstart .+ 4 .* AB + @test (@inferred execute_direct!(plan, 1.0, 0.0)) === plan.Cstorage + @test Array(Cv) ≈ AB +end diff --git a/test/execution/test_halfpack.jl b/test/execution/test_halfpack.jl new file mode 100644 index 0000000..624b98e --- /dev/null +++ b/test/execution/test_halfpack.jl @@ -0,0 +1,495 @@ +# execute_half_packed!: the opt-in "pack B, read a one-tile A in place" path +# (src/execution/halfpack.jl). Two properties are checked: +# +# 1. `_half_pack_a_eligible` admits exactly the plans it should: real +# eltype, SIMDKernel, dense A storage, Qm == mr(kernel), unit-stride A rows. +# 2. `execute_half_packed!` is NEVER wrong. On an eligible plan it must match +# execute! BITWISE (`==`): the per-K-step arithmetic is identical -- same +# `_accumulate_step`, FMA order, K blocking and store, on the same A +# values read from a different address -- and that is what is observed. +# On an ineligible plan it IS execute!, so `==` holds by construction. +# execute_tilewise! (the oracle) and a brute-force reference are +# compared too: `==` with the oracle on eligible plans (observed as well, +# same kernel and K panels), `≈` with the reference. +# +# Every run uses its own freshly built plan and operands, so no path can see +# another's output or buffers. On the half-packed run the plan's packed-A +# buffer is poisoned with NaN first: an eligible run must leave it untouched +# (A really is not packed) and still produce a finite, correct result. + +const execute_half_packed! = QuasiStrided.execute_half_packed! +const _hp_eligible = QuasiStrided._half_pack_a_eligible + +# A DenseMatrix that is not a DenseVector: StridedViews keeps any non-`Array` +# parent as-is (`_normalizeparent`), so a plan over it has `Astorage` of this +# type, which the half-packed path's dense-storage clause must reject. +struct _HPDenseMat{T} <: DenseMatrix{T} + data::Matrix{T} +end +Base.size(a::_HPDenseMat) = size(a.data) +Base.IndexStyle(::Type{<:_HPDenseMat}) = IndexLinear() +Base.getindex(a::_HPDenseMat, i::Int) = a.data[i] +Base.setindex!(a::_HPDenseMat, v, i::Int) = (a.data[i] = v) + +# The default kernel a large-M plan picks for `T` (no extent demotion). +function _hp_default_kernel(::Type{T}) where {T} + plan = plan_contract( + StridedView(zeros(T, 512, 8)), StridedView(zeros(T, 512, 64)), (1, 2), + StridedView(zeros(T, 64, 8)), (2, 3), (1, 3) + ) + return plan.kernel +end + +# Brute-force reference, same as test_direct.jl's `_direct_ref` but kept +# separate so this file does not depend on another test file's helpers. +function _hp_ref(A, indA, B, indB, Cstart, indC, alpha, beta; conjA = false, conjB = false) + dims = Dict{Int, Int}() + for (l, L) in zip(indA, size(A)) + dims[l] = L + end + for (l, L) in zip(indB, size(B)) + dims[l] = L + end + klabels = Tuple(l for l in indA if l in indB && !(l in indC)) + ksizes = Tuple(dims[l] for l in klabels) + T = eltype(Cstart) + out = similar(Cstart, T) + for Ic in CartesianIndices(Cstart) + at = Dict{Int, Int}(zip(indC, Tuple(Ic))) + acc = zero(T) + for Ik in CartesianIndices(ksizes) + for (l, i) in zip(klabels, Tuple(Ik)) + at[l] = i + end + a = A[(at[l] for l in indA)...] + b = B[(at[l] for l in indB)...] + acc += (conjA ? conj(a) : a) * (conjB ? conj(b) : b) + end + out[Ic] = iszero(beta) ? alpha * acc : alpha * acc + beta * Cstart[Ic] + end + return out +end + +# Run execute!, execute_tilewise! and execute_half_packed! on three fresh +# fixtures from `mk()` (identical contents on every call). Returns +# `(C_exec, C_tw, C_half, eligible)`, where `eligible` is +# `_half_pack_a_eligible` of the half-packed run's plan. +function _hp_three_way(mk, alpha, beta; plankw...) + out = Any[] + eligible = false + for run! in (execute!, execute_tilewise!, execute_half_packed!) + Cv, Av, indA, Bv, indB, indC = mk() + plan = plan_contract(Cv, Av, indA, Bv, indB, indC; plankw...) + if run! === execute_half_packed! + eligible = _hp_eligible(plan) + fill!(plan.workspace.packed_a, NaN) + end + ret = run!(plan, alpha, beta) + @test ret === plan.Cstorage + if run! === execute_half_packed! && eligible + # A was read in place, never packed. + @test all(isnan, plan.workspace.packed_a) + end + push!(out, Array(Cv)) + end + return (out..., eligible) +end + +# Column-major matmul fixture C[m,n] = A[m,k] * B[k,n] with a choice of A +# layout. Gaps in every padded parent are NaN, so reading one poisons C. +# :dense plain Matrix (unit-stride rows) +# :ldgap rows 3:Ma+2 of a taller NaN-padded matrix, nonzero offset +# (unit-stride rows) +# :transposed a Ka-by-Ma matrix viewed transposed (row stride Ka) +# :rowgap every other row of a 2Ma-row matrix (row stride 2) +# :wrapped `_HPDenseMat` parent (non-DenseVector storage) +function _hp_mm_maker(::Type{T}, Ma, Ka, Na, seed; alayout = :dense, Cfill = nothing) where {T} + return function () + rng = MersenneTwister(seed) + Amat = randn(rng, T, Ma, Ka) + Bmat = randn(rng, T, Ka, Na) + Cmat = Cfill === nothing ? randn(rng, T, Ma, Na) : fill(convert(T, Cfill), Ma, Na) + Av = if alayout === :dense + StridedView(Amat) + elseif alayout === :ldgap + big = fill(convert(T, NaN), Ma + 5, Ka) + big[3:(Ma + 2), :] .= Amat + StridedView(view(big, 3:(Ma + 2), :)) + elseif alayout === :transposed + permutedims(StridedView(permutedims(Amat, (2, 1))), (2, 1)) + elseif alayout === :rowgap + big = fill(convert(T, NaN), 2Ma, Ka) + big[1:2:(2Ma), :] .= Amat + StridedView(view(big, 1:2:(2Ma), :)) + elseif alayout === :wrapped + StridedView(_HPDenseMat(copy(Amat)), (Ma, Ka), (1, Ma), 0) + else + error("unknown alayout $alayout") + end + @assert Array(Av) == Amat + return (StridedView(Cmat), Av, (1, 2), StridedView(Bmat), (2, 3), (1, 3)) + end +end + +# The (alpha, beta) pairs every eligible shape runs: overwrite, general, +# accumulate, negative scale. +const _HP_AB = ((1.0, 0.0), (2.5, -0.75), (1.0, 1.0), (-0.5, 1.25)) + +@testset "half-pack eligibility: which plans take the in-place A path ($T)" for T in (Float64, Float32) + kern = _hp_default_kernel(T) + MR = mr(kern) + planof(mk; kw...) = (f = mk(); plan_contract(f[1], f[2], f[3], f[4], f[5], f[6]; kw...)) + + # Eligible: M exactly one register tile, unit-stride A rows in dense storage. + for K in (1, 2, 7, 40), N in (1, 5, 13) + p = planof(_hp_mm_maker(T, MR, K, N, 1)) + @test mr(p.kernel) == MR + @test _hp_eligible(p) + @test _hp_eligible(planof(_hp_mm_maker(T, MR, K, N, 1; alayout = :ldgap))) + end + # conjA/conjB on a real eltype: plan_contract folds both to `identity` + # (src/planning/conjugation.jl, `_qs_isconj`), so the plan stays eligible. + for conjA in (false, true), conjB in (false, true) + p = planof(_hp_mm_maker(T, MR, 7, 5, 1); conjA = conjA, conjB = conjB) + @test p.atransform === identity && p.btransform === identity + @test _hp_eligible(p) + end + # ... including explicitly named SIMD kernels of other shapes. + for k in (SIMDKernel(Val(8), Val(6), T), SIMDKernel(Val(2 * lanewidth(kern)), Val(3), T)) + @test _hp_eligible(planof(_hp_mm_maker(T, mr(k), 9, 7, 2); kernel = k)) + @test !_hp_eligible(planof(_hp_mm_maker(T, mr(k) + 1, 9, 7, 2); kernel = k)) + end + + # Ineligible: M one short / one over the tile (whatever kernel the planner + # then picks, Qm != mr), non-unit-stride rows, non-DenseVector storage, + # a ScalarKernel. + @test !_hp_eligible(planof(_hp_mm_maker(T, MR - 1, 7, 5, 3))) + @test !_hp_eligible(planof(_hp_mm_maker(T, MR + 1, 7, 5, 3))) + @test !_hp_eligible(planof(_hp_mm_maker(T, MR, 7, 5, 3; alayout = :transposed))) + @test !_hp_eligible(planof(_hp_mm_maker(T, MR, 7, 5, 3; alayout = :rowgap))) + pw = planof(_hp_mm_maker(T, MR, 7, 5, 3; alayout = :wrapped)) + @test !(pw.Astorage isa DenseVector) + @test !_hp_eligible(pw) + @test !_hp_eligible(planof(_hp_mm_maker(T, 8, 7, 5, 3); kernel = ScalarKernel(Val(8), Val(4), T))) +end + +@testset "half-pack eligibility: complex plans are always excluded ($T)" for T in (ComplexF64, ComplexF32) + for conjA in (false, true), conjB in (false, true), wrap in (false, true) + # Build once to learn the kernel's mr, then again with M == mr so that + # the element type is the only clause that can fail. + p0 = plan_contract( + StridedView(zeros(T, 512, 8)), StridedView(zeros(T, 512, 8)), (1, 2), + StridedView(zeros(T, 8, 8)), (2, 3), (1, 3) + ) + MR = mr(p0.kernel) + Cv, Av, iA, Bv, iB, iC = _hp_mm_maker(T, MR, 6, 5, 4)() + p = plan_contract(Cv, wrap ? conj(Av) : Av, iA, Bv, iB, iC; conjA = conjA, conjB = conjB) + @test axis_length(p.mgroup) == mr(p.kernel) + @test !_hp_eligible(p) + end +end + +@testset "execute_half_packed!: eligible matmul shapes, bitwise vs execute!/oracle ($T)" for T in (Float64, Float32) + MR = mr(_hp_default_kernel(T)) + shapes = ((1, 1), (1, 7), (5, 1), (2, 6), (7, 13), (40, 29), (63, 2), (300, 19)) + for (idx, (K, N)) in enumerate(shapes), alayout in (:dense, :ldgap), (alpha, beta) in _HP_AB + mk = _hp_mm_maker(T, MR, K, N, 100 + idx; alayout = alayout) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta) + @test el + @test C_half == C_exec + @test C_half == C_tw + Cv0, Av0, _, Bv0, _, _ = mk() + @test C_half ≈ _hp_ref(Array(Av0), (1, 2), Array(Bv0), (2, 3), Array(Cv0), (1, 3), alpha, beta) + end +end + +@testset "execute_half_packed!: several jc/pc blocks and explicit kernels ($T)" for T in (Float64, Float32) + kern = _hp_default_kernel(T) + kernels = (kern, SIMDKernel(Val(8), Val(6), T), SIMDKernel(Val(2 * lanewidth(kern)), Val(3), T)) + # (kc, nc): K and N both split into several blocks, with tails. + for k in kernels, (kc, nc) in ((3, 6), (5, 7), (64, 1), (1, 100)) + for (K, N) in ((11, 17), (4, 25), (1, 3)) + mk = _hp_mm_maker(T, mr(k), K, N, 7 * K + N) + C_exec, C_tw, C_half, el = _hp_three_way(mk, 1.5, 0.5; kernel = k, kc = kc, nc = nc) + @test el + @test C_half == C_exec + @test C_half == C_tw + Cv0, Av0, _, Bv0, _, _ = mk() + @test C_half ≈ _hp_ref(Array(Av0), (1, 2), Array(Bv0), (2, 3), Array(Cv0), (1, 3), 1.5, 0.5) + end + end +end + +@testset "execute_half_packed!: K = 0, alpha = 0, beta = 0 never reads what it must not ($T)" for T in (Float64, Float32) + MR = mr(_hp_default_kernel(T)) + N = 9 + # K = 0: beta applied once, A/B never read. + for beta in (0.5, 1.0, 0.0, -2.0) + Cstart = randn(MersenneTwister(31), T, MR, N) + mk = () -> ( + StridedView(copy(Cstart)), StridedView(fill(convert(T, NaN), MR, 0)), (1, 2), + StridedView(fill(convert(T, Inf), 0, N)), (2, 3), (1, 3), + ) + C_exec, C_tw, C_half, _ = _hp_three_way(mk, 1.0, beta) + @test C_half == C_exec + @test C_half == C_tw + @test C_half ≈ beta .* Cstart + end + # alpha = 0: beta applied once, A/B (NaN/Inf) never read. + for beta in (1.25, 1.0, 0.0) + Cstart = randn(MersenneTwister(32), T, MR, N) + mk = () -> ( + StridedView(copy(Cstart)), StridedView(fill(convert(T, NaN), MR, 6)), (1, 2), + StridedView(fill(convert(T, Inf), 6, N)), (2, 3), (1, 3), + ) + C_exec, C_tw, C_half, el = _hp_three_way(mk, 0.0, beta) + @test el + @test C_half == C_exec + @test all(isfinite, C_half) + @test C_half ≈ beta .* Cstart + end + # beta = 0 into a NaN-poisoned C: old C never read. + for (K, N2) in ((1, 1), (5, 4), (40, 13)) + mk = _hp_mm_maker(T, MR, K, N2, 900 + K; Cfill = NaN) + C_exec, C_tw, C_half, el = _hp_three_way(mk, 1.5, 0.0) + @test el + @test all(isfinite, C_half) + @test C_half == C_exec + @test C_half == C_tw + end + # Empty N: a no-op. + Cmat = fill(convert(T, NaN), MR, 0) + plan = plan_contract( + StridedView(Cmat), StridedView(randn(T, MR, 4)), (1, 2), + StridedView(randn(T, 4, 0)), (2, 3), (1, 3) + ) + @test execute_half_packed!(plan, 1.0, 2.0) === plan.Cstorage +end + +@testset "execute_half_packed!: scattered K, negative-stride B, offset C, non-ramp M ($T)" for T in (Float64, Float32) + kern = _hp_default_kernel(T) + MR = mr(kern) + + # K composite (k1, k2) that does not fold in A (a gapped slice), so K is a + # scattered axis on the in-place A read; B reversed along k1 (negative + # stride); C a slice with a nonzero offset. Every gap in A's parent is NaN. + for N in (1, 6, 23), (kc, nc) in ((64, 64), (3, 6)) + mk = function () + # `local`: these names are also locals of the enclosing loop, and + # a closure would otherwise assign the outer ones. + local Abig, Av, Bv, Cbig, Cv + rng = MersenneTwister(50 + N) + Abig = fill(convert(T, NaN), MR + 3, 5, 9) + Abig[2:(MR + 1), 1:4, 2:2:9] .= randn(rng, T, MR, 4, 4) + Av = StridedView(view(Abig, 2:(MR + 1), 1:4, 2:2:9)) + Bv = StridedView(view(randn(rng, T, 4, 4, N), 4:-1:1, :, :)) + Cbig = randn(rng, T, MR + 2, N + 3) + Cv = StridedView(view(Cbig, 2:(MR + 1), 2:(N + 1))) + return (Cv, Av, (1, 2, 3), Bv, (2, 3, 4), (1, 4)) + end + Cv, Av, iA, Bv, iB, iC = mk() + p = plan_contract(Cv, Av, iA, Bv, iB, iC; kernel = kern) + @test !QuasiStrided.affine_ramp(p.kgroup)[1] + @test any(<(0), strides(Bv)) && offset(Cv) != 0 + for (alpha, beta) in ((1.0, 0.0), (-1.5, 0.5)) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta; kernel = kern, kc = kc, nc = nc) + @test el + @test C_half == C_exec + @test C_half == C_tw + @test all(isfinite, C_half) + @test C_half ≈ _hp_ref(Array(Av), iA, Array(Bv), iB, Array(Cv), iC, alpha, beta) + end + end + + # M composite (a1, a2) contiguous in A but gapped in C (C[a1, n, a2]): + # not an affine ramp, so the M descriptors are materialized and C's rows + # are a scattered axis, while A's rows are still one unit-stride run. + for N in (2, 7) + mk = function () + local Av, Bv, Cv + rng = MersenneTwister(60 + N) + Av = StridedView(randn(rng, T, MR ÷ 2, 2, 5)) + Bv = StridedView(randn(rng, T, 5, N)) + Cv = StridedView(randn(rng, T, MR ÷ 2, N, 2)) + return (Cv, Av, (1, 2, 3), Bv, (3, 4), (1, 4, 2)) + end + Cv, Av, iA, Bv, iB, iC = mk() + p = plan_contract(Cv, Av, iA, Bv, iB, iC; kernel = kern) + @test !QuasiStrided.affine_ramp(p.mgroup)[1] + C_exec, C_tw, C_half, el = _hp_three_way(mk, 0.75, -1.0; kernel = kern) + @test el + @test C_half == C_exec + @test C_half == C_tw + @test C_half ≈ _hp_ref(Array(Av), iA, Array(Bv), iB, Array(Cv), iC, 0.75, -1.0) + end + + # helpers.jl's `scattered_fixture` (benchmark/harness.jl's `build_scattered` + # shape: permuted A with a zero-stride b axis, negative-stride B, sliced + # C) with M forced to one register tile. b_n = 1: M is `a` alone, unit + # stride in A -> eligible. b_n = 2: M = (a, b) with b's A stride 0, so A's + # rows repeat -> not affine -> ineligible, execute! fallback. + sentinel = convert(T, 7) + for (a_n, b_n, want) in ((MR, 1, true), (MR ÷ 2, 2, false)), (k_n, n_n) in ((32, 32), (3, 5), (1, 1)) + Cstart = randn(MersenneTwister(40 + a_n + k_n), T, a_n, n_n, b_n) + mk = function () + local Cv, Av, indA, Bv, indB, indC + Cv, Av, indA, Bv, indB, indC = scattered_fixture(T, a_n, k_n, b_n, n_n) + fill!(parent(Cv), sentinel) + copyto!(Cv, Cstart) + return (Cv, Av, indA, Bv, indB, indC) + end + for (alpha, beta) in ((1.0, 0.0), (-1.5, 0.5)) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta) + @test el == want + @test C_half == C_exec + want && @test C_half == C_tw + @test C_half ≈ C_tw + Cv, Av, indA, Bv, indB, indC = mk() + @test C_half ≈ _hp_ref(Array(Av), indA, Array(Bv), indB, Cstart, indC, alpha, beta) + # Nothing outside C's slice is written. + plan = plan_contract(Cv, Av, indA, Bv, indB, indC) + execute_half_packed!(plan, alpha, beta) + @test Array(Cv) == C_half + fill!(Cv, sentinel) + @test all(==(sentinel), parent(Cv)) + end + end +end + +@testset "execute_half_packed!: operand-role swap ($T)" for T in (Float64, Float32) + # C[n,m] column-major with n of extent MR: plan_contract swaps roles, so + # the plan's M (and its in-place "A") is the user's B. It is eligible iff + # the user's B is unit-stride along n. + MR = mr(_hp_default_kernel(T)) + Ma, Ka = 5, 7 + for btrans in (true, false) + mk = function () + local Amat, Bmat, Cmat, Bv + rng = MersenneTwister(606) + Amat = randn(rng, T, Ma, Ka) + Bmat = randn(rng, T, Ka, MR) + Cmat = randn(rng, T, MR, Ma) + Bv = btrans ? permutedims(StridedView(permutedims(Bmat, (2, 1))), (2, 1)) : StridedView(Bmat) + return (StridedView(Cmat), StridedView(Amat), (1, 2), Bv, (2, 3), (3, 1)) + end + Cv, Av, iA, Bv, iB, iC = mk() + plan = plan_contract(Cv, Av, iA, Bv, iB, iC) + @test plan.Astorage === parent(Bv) + for (alpha, beta) in ((1.0, 0.0), (0.5, -1.0)) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta) + @test el == btrans + @test C_half == C_exec + @test C_half ≈ C_tw + @test C_half ≈ _hp_ref(Array(Av), iA, Array(Bv), iB, Array(Cv), iC, alpha, beta) + end + end +end + +@testset "execute_half_packed!: ineligible plans fall back to execute! exactly ($T)" for T in (Float64, Float32) + MR = mr(_hp_default_kernel(T)) + cases = ( + (MR - 1, 7, 5, :dense, (;)), + (MR + 1, 7, 5, :dense, (;)), + (2MR + 3, 9, 4, :dense, (;)), + (MR, 7, 5, :transposed, (;)), + (MR, 7, 5, :rowgap, (;)), + (MR, 7, 5, :wrapped, (;)), + (8, 7, 5, :dense, (kernel = ScalarKernel(Val(8), Val(4), T),)), + (8, 7, 5, :dense, (kernel = ScalarKernel(Val(8), Val(4), T), kc = 3, nc = 4)), + ) + for (idx, (M, K, N, alayout, kw)) in enumerate(cases), (alpha, beta) in _HP_AB + mk = _hp_mm_maker(T, M, K, N, 700 + idx; alayout = alayout) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta; kw...) + @test !el + @test C_half == C_exec + @test C_half ≈ C_tw + Cv0, Av0, _, Bv0, _, _ = mk() + @test C_half ≈ _hp_ref(Array(Av0), (1, 2), Array(Bv0), (2, 3), Array(Cv0), (1, 3), alpha, beta) + end +end + +@testset "execute_half_packed!: complex plans fall back to execute! exactly ($T)" for T in (ComplexF64, ComplexF32) + p0 = plan_contract( + StridedView(zeros(T, 512, 8)), StridedView(zeros(T, 512, 8)), (1, 2), + StridedView(zeros(T, 8, 8)), (2, 3), (1, 3) + ) + MR = mr(p0.kernel) + alpha, beta = (0.75 - 1.0im), (0.5 + 0.5im) + for (M, K, N) in ((MR, 1, 1), (MR, 4, 3), (MR, 9, 13), (MR - 1, 5, 2), (3, 4, 2)), + conjA in (false, true), conjB in (false, true) + mk = _hp_mm_maker(T, M, K, N, 3000 + M + K) + C_exec, C_tw, C_half, el = _hp_three_way(mk, alpha, beta; conjA = conjA, conjB = conjB) + @test !el + @test C_half == C_exec + @test C_half ≈ C_tw + Cv, Av, _, Bv, _, _ = mk() + @test C_half ≈ _hp_ref( + Array(Av), (1, 2), Array(Bv), (2, 3), Array(Cv), (1, 3), alpha, beta; + conjA = conjA, conjB = conjB + ) + # A conj-wrapped view composes with the flag by XOR. + Cvv, Avv, iA, Bvv, iB, iC = mk() + planv = plan_contract(Cvv, conj(Avv), iA, Bvv, iB, iC; conjA = !conjA, conjB = conjB) + @test !_hp_eligible(planv) + execute_half_packed!(planv, alpha, beta) + @test Array(Cvv) == C_half + end +end + +@testset "execute_half_packed!: out-of-bounds operands are rejected like execute! ($T)" for T in (Float64, Float32) + # StridedView does not validate its extent against the parent, and + # plan_contract does not either; the drivers' hoisted span checks do. The + # in-place A loads are raw-pointer reads, so the A check is load-bearing. + MR = mr(_hp_default_kernel(T)) + K, N = 6, 5 + mkA = () -> ( + StridedView(zeros(T, MR, N)), StridedView(randn(T, MR * K), (MR, K), (1, MR), 3), (1, 2), + StridedView(randn(T, K, N)), (2, 3), (1, 3), + ) + mkB = () -> ( + StridedView(zeros(T, MR, N)), StridedView(randn(T, MR, K)), (1, 2), + StridedView(randn(T, K * N), (K, N), (1, K), 1), (2, 3), (1, 3), + ) + mkC = () -> ( + StridedView(zeros(T, MR * N), (MR, N), (1, MR), 2), StridedView(randn(T, MR, K)), (1, 2), + StridedView(randn(T, K, N)), (2, 3), (1, 3), + ) + for mk in (mkA, mkB, mkC), run! in (execute!, execute_half_packed!) + Cv, Av, iA, Bv, iB, iC = mk() + plan = plan_contract(Cv, Av, iA, Bv, iB, iC) + @test _hp_eligible(plan) + @test_throws BoundsError run!(plan, 1.0, 0.0) + end +end + +@testset "execute_half_packed!: allocation-free, reusable, inferred ($T)" for T in (Float64, Float32) + kern = _hp_default_kernel(T) + MR = mr(kern) + # Affine K and scattered K (the `PtrScatterAxis` A view). + Cmat = zeros(T, MR, 37) + plan = _mm_plan(Cmat, randn(T, MR, 50), randn(T, 50, 37); kc = 16, nc = 12) + @test _hp_eligible(plan) + @test _steady_allocs!(execute_half_packed!, plan, Cmat) == 0 + + Abig = randn(T, MR + 3, 5, 9) + Csc = zeros(T, MR, 11) + plan_sc = plan_contract( + StridedView(Csc), StridedView(view(Abig, 2:(MR + 1), 1:4, 2:2:9)), (1, 2, 3), + StridedView(randn(T, 4, 4, 11)), (2, 3, 4), (1, 4); kernel = kern, kc = 5 + ) + @test _hp_eligible(plan_sc) + @test !QuasiStrided.affine_ramp(plan_sc.kgroup)[1] + @test _steady_allocs!(execute_half_packed!, plan_sc, Csc) == 0 + + # A plan is reusable: a second call accumulates onto the first's output. + Amat, Bmat = randn(T, MR, 9), randn(T, 9, 8) + Cstart = randn(T, MR, 8) + C = copy(Cstart) + p = _mm_plan(C, Amat, Bmat) + @test _hp_eligible(p) + execute_half_packed!(p, 2, 1) + execute_half_packed!(p, 2, 1) + @test C ≈ Cstart .+ 4 .* (Amat * Bmat) + @test (@inferred execute_half_packed!(p, 1.0, 0.0)) === p.Cstorage + @test C ≈ Amat * Bmat +end diff --git a/test/helpers.jl b/test/helpers.jl index ca8114b..607736a 100644 --- a/test/helpers.jl +++ b/test/helpers.jl @@ -4,8 +4,9 @@ using StridedViews: StridedView, offset # plan_contract/execute!/ContractPlan aren't exported (only contract! is); -# execute_tilewise! is never exported at all (it's an internal oracle). -# `test/runtests.jl` deliberately leaves these four out of its own +# execute_tilewise! is never exported at all (it's an internal oracle), nor is +# the opt-in execute_direct!. +# `test/runtests.jl` deliberately leaves these five out of its own # name-restoring `using QuasiStrided: ...` block, because a `const` may not # shadow an imported binding; the workspace API is reached as # `QuasiStrided.` for the same reason. @@ -13,6 +14,7 @@ const plan_contract = QuasiStrided.plan_contract const execute! = QuasiStrided.execute! const ContractPlan = QuasiStrided.ContractPlan const execute_tilewise! = QuasiStrided.execute_tilewise! +const execute_direct! = QuasiStrided.execute_direct! # TO names are always qualified: a bare `using TensorOperations` collides with # QuasiStrided's `scalartype`. diff --git a/test/microkernels/test_fmaddsub_kernel.jl b/test/microkernels/test_fmaddsub_kernel.jl index badc682..df44b0b 100644 --- a/test/microkernels/test_fmaddsub_kernel.jl +++ b/test/microkernels/test_fmaddsub_kernel.jl @@ -541,8 +541,18 @@ const _QSF = QuasiStrided # shape, read off the hot loop of `accumulate` with the driver's # `PackedPanel` argument types (benchmark/probes/fmaddsub_codegen.jl # has the same count for every shape, including under `-C znver2`). + # + # Also skipped on CI (GitHub Actions sets CI=true): these assert an + # EXACT instruction mix (down to spill counts) in the compiled hot + # loop, which depends on the host's real, physically-exposed CPU + # features -- not something a generic hosted CI runner reliably + # matches (observed: mismatched counts on both ubuntu-latest and + # macos-latest runners, presumably a virtualized/hypervisor-exposed + # feature set that doesn't line up with this test's tuned + # expectations). This is exactly the kind of measurement this + # project otherwise takes on dedicated Slurm hardware, not shared CI. isa = target_profile().isa - has_fma = Sys.ARCH === :x86_64 && isa in (:avx2, :avx512) + has_fma = Sys.ARCH === :x86_64 && isa in (:avx2, :avx512) && get(ENV, "CI", "false") != "true" function hot_loop(asm) lines = split(asm, '\n') labels = Dict{String, Int}() diff --git a/test/runtests.jl b/test/runtests.jl index d2c5fe0..cafc495 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,6 +49,8 @@ using QuasiStrided: AxisGroup, axis_length, offsets, fill_offsets!, BlockDescrip include("execution/test_workspace.jl") include("execution/test_scalar_vs_simd.jl") include("execution/test_macro_blocking.jl") + include("execution/test_direct.jl") + include("execution/test_halfpack.jl") # Last among the functional tests: it does a bare `using TensorOperations`, # which exports its own `scalartype` and would otherwise conflict with the