Performance and simplification: per-call floor, scrambled layouts, AVX-512 GEMM peak, degenerate shapes, compile latency

README.md

@@ -5,15 +5,12 @@
 
 A native-Julia dense tensor contraction engine for `StridedView`s
 (`StridedViews.jl`): grouped-axis / block-scatter indexing, packing, and
-fixed-shape microkernels, in the style of block-scatter-matrix tensor
+fixed-shape SIMD microkernels, in the style of block-scatter-matrix tensor
 contraction (BSMTC, Matthews arXiv:1607.00291). It ships as a
 [TensorOperations.jl](https://github.com/QuantumKitHub/TensorOperations.jl)
 backend, `QuasiStridedBackend`.
 
-**Status**: not registered. Includes a BLIS five-loop (`NC`/`KC`/`MC`)
-macro-blocking driver with packed-panel reuse, on top of a scalar reference
-and a SIMD kernel (the engine-wide default), plus planar and 1m complex
-microkernels for `ComplexF32`/`ComplexF64`. See "Status" below.
+**Status**: not registered; see [Status](#status).
 
 ## Install
 
@@ -24,9 +21,7 @@ Pkg.add(url="https://github.com/lkdvos/QuasiStrided.jl")
 
 ## Usage
 
-The primary, intended way to use this package is through TensorOperations'
-own `@tensor`/`ncon` macros, opting into QuasiStrided's engine with
-`backend=QuasiStridedBackend()`:
+Opt into the engine through TensorOperations' `@tensor`/`ncon`:
 
 ```julia
 using TensorOperations, QuasiStrided
@@ -37,165 +32,109 @@ A, B = randn(3, 5, 2), randn(5, 4)
 @tensor backend = QuasiStridedBackend() C[a, n, b] := A[a, k, b] * B[k, n]
 
 # equivalently, via ncon
-C2 = ncon((A, B), ((-1, 1, -3), (1, -2)); backend = QuasiStridedBackend())
+C2 = ncon([A, B], [[-1, 1, -3], [1, -2]]; backend = QuasiStridedBackend())
 ```
 
-Loading QuasiStrided never changes the behavior of plain `@tensor` code —
-`QuasiStridedBackend` is *not* hooked into `TensorOperations.select_backend`,
-so it only ever runs when named explicitly. See "Not implemented" below for
-what it accepts and rejects.
+`QuasiStridedBackend` is not registered with `TensorOperations.select_backend`,
+so loading QuasiStrided never changes the behaviour of `@tensor` code that does
+not name it.
 
-### Engine interface (direct/low-level usage)
+### Engine interface
 
-`QuasiStridedBackend` is a thin adapter over QuasiStrided's own contraction
-engine, which remains available directly (and is what the TensorOperations
-backend calls internally):
+The engine behind the backend can also be called directly:
 
 ```julia
 using QuasiStrided, StridedViews
+using QuasiStrided: contract!, plan_contract, execute!, SIMDKernel
 
 # C[a,n,b] = alpha * sum_k A[a,k,b] * B[k,n] + beta * C[a,n,b]
 A, B, C = randn(3, 5, 2), randn(5, 4), zeros(3, 4, 2)
 
-# indA/indB/indC label each axis; a label shared by A and B but absent from
-# indC is contracted, one present in indC is kept (M if from A, N if from B).
+# A label shared by A and B but absent from indC is contracted.
 contract!(StridedView(C), 1.0,
           StridedView(A), (1, 2, 3),   # a, k, b
           StridedView(B), (2, 4),      # k, n
           0.0,
           (1, 4, 3))                   # a, n, b
 
-# Reuse a plan/workspace across repeated calls of the same shape, and pick a
-# kernel explicitly (default: SIMDKernel):
+# Reuse a plan and its workspace across calls, optionally naming the kernel:
 plan = plan_contract(StridedView(C), StridedView(A), (1, 2, 3),
-                      StridedView(B), (2, 4), (1, 4, 3);
-                      kernel=SIMDKernel(Val(8), Val(6), Float64))
+                     StridedView(B), (2, 4), (1, 4, 3);
+                     kernel = SIMDKernel(Val(8), Val(6), Float64))
 execute!(plan, 1.0, 0.0)
 ```
 
-`plan_contract` also accepts `workspace=`/`allocator=` keywords to control
-buffer reuse and, e.g., arena/Bumper-backed allocation instead of the
-default GC-owned, task-local pool; see the docstring for details.
+`plan_contract` also takes `workspace=`/`allocator=` keywords (e.g. for
+Bumper-backed buffers); see its docstring.
 
 ## API
 
-QuasiStrided exports exactly one name:
+The only export is `QuasiStridedBackend`. The following names are `public`
+(Julia >= 1.11) but unexported, to avoid collisions with TensorOperations
+(both define e.g. `scalartype`):
 
-| Name | Kind |
+| Kind | Names |
 | --- | --- |
-| `QuasiStridedBackend` | Exported |
+| Contraction | `contract!`, `plan_contract`, `execute!`, `ContractPlan`, `ContractWorkspace` |
+| Blocking | `Blocking`, `default_blocking` |
+| Kernels | `ScalarKernel`, `SIMDKernel`, `PlanarKernel`, `OneMKernel`, `FMAddSubKernel` |
+| Hardware | `target_profile`, `cache_topology`, `TargetProfile`, `CacheLevel` |
 
-Everything else the engine needs to be used directly is `public` but
-unexported (Julia >= 1.11; visible via `QuasiStrided.<name>` on 1.10):
-
-| Public, unexported names |
-| --- |
-| `contract!`, `plan_contract`, `execute!`, `ContractPlan`, `ContractWorkspace`, `Blocking`, `default_blocking`, `ScalarKernel`, `SIMDKernel`, `PlanarKernel`, `OneMKernel`, `target_profile`, `cache_topology`, `TargetProfile`, `CacheLevel` |
-
-Everything not listed above (indexing, packing, tile/kernel internals such
-as `AxisGroup`, `KernelDescriptor`, `pack_a!`/`pack_b!`, `accumulate`, etc.)
-is internal implementation detail: it carries a docstring and is freely used
-within the package and its tests, but no semver promise and no expectation
-of stability across releases. The split has three tiers (1 exported + 15
-public + internal names); keeping most names unexported avoids export
-collisions with TensorOperations (e.g. both packages have a `scalartype`).
+Everything else is internal.
 
 ## Code map
 
-`src/` follows the algorithm. `plan_contract` runs the planning stages once;
-`execute!` runs a BLIS five-loop nest that packs slivers of `A` and `B` into
-contiguous panels and drives a microkernel over register tiles of `C`.
+`plan_contract` resolves labels into M/N/K axis groups and picks a kernel,
+blocking and workspace; `execute!` runs a BLIS five-loop nest that packs `A`
+and `B` into panels and drives a microkernel over register tiles of `C`.
 
-| Folder | Stage | Contents |
-| --- | --- | --- |
-| `src/hardware/` | detection | ISA, vector width, register count and cache topology, detected once per process (`target_profile()`) |
-| `src/layout/` | addressing | `AxisGroup` (a grouped tensor axis as zero-based offsets), block descriptors, and the affine/scattered tile axes and `QSTile` the packers and kernels read and write through |
-| `src/packing/` | packing | the packed-panel formats and the `Descriptor` that fixes them per kernel (`format.jl`), borrowed-pointer panels, and the `pack_a!`/`pack_b!` packers with their contiguous fast paths |
-| `src/microkernels/` | microkernels | the shared kernel interface (`execute_tile!` = `zero_accumulator` + `accumulate` + `store_tile!`) and the scalar, SIMD, planar-complex and 1m-complex kernels |
-| `src/planning/` | planning | label classification and ordering, conjugation, kernel selection from the detected hardware, cache blocking, and `plan_contract`/`ContractPlan` |
-| `src/execution/` | execution | the reusable `ContractWorkspace`, the five-loop nest behind `execute!`/`contract!`, its macro-kernel helpers, and `execute_tilewise!`, the tile-by-tile correctness oracle |
-| `src/integrations/` | adapters | `QuasiStridedBackend` for TensorOperations.jl |
+| Folder | Contents |
+| --- | --- |
+| `src/hardware/` | ISA, vector width, register count and cache topology, detected once per process |
+| `src/layout/` | `AxisGroup` (a grouped axis as zero-based offsets), block descriptors, tile axes |
+| `src/packing/` | packed-panel formats and the `pack_a!`/`pack_b!` packers with contiguous fast paths |
+| `src/microkernels/` | the kernel interface and the scalar, SIMD, planar, 1m and fmaddsub kernels |
+| `src/planning/` | label classification, conjugation, kernel selection, cache blocking, `plan_contract` |
+| `src/execution/` | the workspace, the five-loop nest, the specialised paths and the tile-by-tile oracle |
+| `src/integrations/` | `QuasiStridedBackend` for TensorOperations.jl |
 
-`test/` mirrors this layout; `test/runtests.jl` includes every file into one
-scope, in stage order.
+`test/` mirrors this layout.
 
 ## Status
 
-The package implements:
-
-- Grouped-axis (block-scatter) indexing: each tensor axis group is a list of
-  zero-based offsets, checked for overflow, over a `StridedView`.
-- Packing of `A` and `B` slivers into kernel-specific panel formats, fixed per
-  kernel by one descriptor, with contiguous fast paths.
-- Four microkernels: `ScalarKernel` (the reference), `SIMDKernel` (the default
-  for real element types), `PlanarKernel` (split-complex, the default for
-  complex) and `OneMKernel` (the 1m induced method, reusing the real kernel
-  over `2*kc` steps).
-- A BLIS five-loop (`NC`/`KC`/`MC`) macro-blocking `execute!` that reuses a
-  packed B panel per `(jc,pc)` block and a packed A panel per `(jc,pc,ic)`
-  block. `Blocking`/`default_blocking` expose measured `mc`/`kc`/`nc`
-  defaults.
-- A register shape derived from the detected hardware (`target_profile()`:
-  ISA, vector width, register count), with per-ISA overrides and a
-  register-budget fit for ISAs without one. When `M` is too small to fill one
-  register tile, the kernel is demoted to the conservative fitted shape. For
-  real element types and a shallow contracted extent, a kernel whose `MR` does
-  not divide `C`'s leading unit-stride run is demoted to a menu shape that
-  does, so the vectorized store applies.
-- `ComplexF32`/`ComplexF64` support; `conjA`/`conjB` and each operand's
-  `StridedView.op` are folded into packing. A conjugated output is rejected.
-- Zero steady-state allocation for `execute!` on Julia >= 1.11. On Julia 1.10
-  the older compiler does not keep `SIMDKernel`'s accumulator in registers, so
-  it allocates per call; results are still correct.
-- A TensorOperations.jl backend, `QuasiStridedBackend`, with task-local
-  workspace pooling for the default allocator and call-scoped buffers for an
-  explicit `allocator=` (e.g. `Bumper`). `plan_contract`'s per-call
-  bookkeeping still allocates a few KB.
-- `execute_tilewise!`, a tile-by-tile driver used in the tests as the
-  correctness oracle for `execute!`.
-
-**Not implemented** (deliberately):
-
-- **Contraction only, engine-wise.** `QuasiStridedBackend` implements
-  `TensorOperations.tensorcontract!` itself; there is no QuasiStrided analog
-  of a standalone add/permute step or of diagonal/trace support at all in
-  the engine. `tensoradd!`/`tensortrace!` fall back to `StridedNative()` so
-  a `@tensor` network mixing a contraction with an add/permute or trace step
-  can still run wholesale under `backend=QuasiStridedBackend()` -- a timing
-  taken on those two specific operations under this backend measures
-  `StridedNative`, not this engine.
-- **One shared element type, out of four.** `tensorcontract!` accepts inputs
-  that share a single element type out of
-  `Float32`/`Float64`/`ComplexF32`/`ComplexF64` and are all strided; anything
-  else (mixed eltypes, mixed real-and-complex operands, other eltypes, a
-  non-strided array, an aliased output, a conjugated *output* view) raises
-  `ArgumentError` rather than silently falling back to another backend, so a
-  timing taken with `backend=QuasiStridedBackend()` always measures this
-  engine.
-- **Not the default backend.** QuasiStrided defines no
-  `TensorOperations.select_backend` method, so merely loading it never
-  changes the behavior, performance, or allocation profile of existing
-  `@tensor`/`ncon` code that doesn't pass `backend=QuasiStridedBackend()`
-  explicitly. This is deliberate: QuasiStrided is expected to lose to
-  `StridedBLAS()` on most shapes, so silently capturing the default dispatch
-  would be a regression for existing users.
-- einsum string parsing; batch axes, diagonals, or isolated reductions; a
-  separate beta-addend tensor; autotuning across shapes; runtime cache
-  probing (no benchmark ever runs at load or first call; `mc`/`kc`/`nc` are
-  picked by an analytical model of the sysfs/sysctl-detected cache geometry,
-  with fixed constants only where that is undetected -- see `default_blocking`); threading; GPU execution; K padding.
-- On the complex side specifically: the **3m** (Karatsuba) method; **mixed
-  real/complex operands** (promotion belongs in TensorOperations'
-  `promote_contract` layer, not here); writing into a **conjugated output**
-  view; **measured** complex register shapes for AVX2 or NEON (complex
-  *works* there — the engine selects the largest shape that fits the detected
-  register file — but those shapes are fitted, not swept, and no throughput
-  claim is made for them); and the unit-stride plane-to-interleave store, which
-  is measurement-gated and not yet justified.
-
-`benchmark/bench_tensoroperations.jl` is a head-to-head comparison against
-`StridedNative`/`StridedBLAS` on identical shapes. `StridedBLAS()` is faster
-than `QuasiStridedBackend()` on the measured shapes, which is why the backend
-is not hooked into `select_backend`; `QuasiStridedBackend()` in turn beats
-`StridedNative()`. `benchmark/bench_to_suite.jl` runs the same comparison
-over TensorOperations.jl's upstream benchmark suite.
+Implemented:
+
+- Grouped-axis (block-scatter) indexing over `StridedView`s, with packing into
+  kernel-specific panel formats and contiguous fast paths.
+- Microkernels: `ScalarKernel` (reference), `SIMDKernel` (default for real
+  eltypes), `PlanarKernel` (split-complex, default for complex eltypes),
+  `OneMKernel` (the 1m method) and `FMAddSubKernel` (interleaved complex, used
+  for small-M complex contractions on AVX-512). Any of them can be named via
+  `plan_contract(...; kernel = ...)`.
+- Register shapes derived from the detected ISA, demoted when `M` is too short
+  for a full tile or when a smaller tile keeps stores into `C` unit-stride.
+- A BLIS five-loop (`NC`/`KC`/`MC`) nest reusing packed panels, with blocking
+  from an analytical model of the detected cache sizes (fixed constants when
+  they are undetected); plus automatically selected paths for dot-product,
+  outer-product and read-B-in-place shapes.
+- `Float32`/`Float64`/`ComplexF32`/`ComplexF64`; `conjA`/`conjB` and each
+  operand's `StridedView.op` are applied during packing.
+- Zero steady-state allocation on Julia >= 1.11 for `execute!` on a reused
+  plan and for `tensorcontract!` through the backend (which pools workspaces
+  per task). On Julia 1.10 `SIMDKernel`'s accumulator is not kept in
+  registers, so calls allocate; results are unaffected.
+
+Not implemented:
+
+- Anything but contraction: `tensoradd!`/`tensortrace!` under this backend are
+  forwarded to `StridedNative()`, so mixed `@tensor` networks still run.
+- Mixed eltypes (including real with complex; promotion is TensorOperations'
+  `promote_contract` job), other eltypes, non-strided operands, an output
+  aliased with an input, and a conjugated output view: `tensorcontract!` throws
+  an `ArgumentError` instead of falling back to another backend.
+- Being the default backend: `StridedBLAS()` is faster on most shapes.
+- Batch axes, traces/diagonals, threading, GPU execution, the 3m complex
+  method, autotuning.
+
+`benchmark/bench_to_suite.jl` compares `QuasiStridedBackend()` with
+`StridedBLAS()` over TensorOperations.jl's upstream benchmark suite.

benchmark/Project.toml

@@ -3,26 +3,21 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
-ProfileCanvas = "efd6af41-a80b-495e-886c-e51b0c7d77a3"
 QuasiStrided = "8a7b1e6a-8f7e-4c1a-9d5a-3b7f0c9e2d41"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StridedViews = "4db3bf67-4bd7-4b4e-b153-31dc3fb37143"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 TensorOperationsBenchmarks = "d983fc97-4e87-46ba-abd9-b4864e69d4dd"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-# Pinned to PR #303 (QuantumKitHub/TensorOperations.jl#303, branch "benchmark"),
-# commit d1ace864bb2d5a8e5db6d1d6d9d51c14ba19369e, because the PR is unmerged.
-# Follow-up: repoint to the registered TensorOperations.jl release once merged.
+# TensorOperationsBenchmarks is unregistered (QuantumKitHub/TensorOperations.jl#303);
+# repoint to a registered release once it exists.
 [sources]
 QuasiStrided = {path = ".."}
 TensorOperationsBenchmarks = {rev = "d1ace864bb2d5a8e5db6d1d6d9d51c14ba19369e", subdir = "benchmark", url = "https://github.com/QuantumKitHub/TensorOperations.jl"}
 
 [compat]
 CairoMakie = "0.15"
-ProfileCanvas = "0.1"
 StridedViews = "0.5"
 TensorOperations = "5.8"
 julia = "1.11"