Complete the 0.10 release notes and update the NUMA example

docs/src/examples/numa_aware.md

@@ -1,6 +1,6 @@
 # NUMA-aware SAXPY
 
-This example demonstrates how to define and run a SAXPY kernel (single-precision `Y[i] = a * X[i] + Y[i]`) such that it runs efficiently on a system with multiple memory domains ([NUMA](https://en.wikipedia.org/wiki/Non-uniform_memory_access)) using multithreading. (You likely will need to fine-tune the value of `N` on your system of interest if you care about the particular measurement.)
+This example demonstrates how to define and run a SAXPY kernel (single-precision `Y[i] = a * X[i] + Y[i]`) such that it runs efficiently on a system with multiple memory domains ([NUMA](https://en.wikipedia.org/wiki/Non-uniform_memory_access)) with the multithreaded `CPU` backend. (You likely will need to fine-tune the value of `N` on your system of interest if you care about the particular measurement.)
 
 ````@eval
 using Markdown
@@ -15,30 +15,29 @@ $(read(path, String))
 
 **Important remarks:**
 
-1) Pin your threads systematically to the available physical (or virtual) CPU-cores. [ThreadPinning.jl](https://github.com/carstenbauer/ThreadPinning.jl) is your friend.
-2) Opt-out of Julia's dynamic task scheduling (especially task migration) by using `CPU(; static=true)` instead of `CPU()`.
-3) Initialize your data in parallel(!). It is of utmost importance to use a parallel access pattern for initialization that is as similar as possible as the access pattern of your computational kernel. The reason for this is ["NUMA first-touch policy"](https://queue.acm.org/detail.cfm?id=2513149#:~:text=This%20is%20called%20the%20first,policy%20associated%20with%20a%20task.). `KernelAbstractions.zeros(backend, dtype, N)` is your friend.
+The `CPU` backend runs kernels with [PoCL](https://portablecl.org) on PoCL's own threads, not on Julia's. That determines how to apply the usual advice for NUMA systems:
+
+1) Pin the threads that run the kernels. Set PoCL's `POCL_AFFINITY=1` environment variable before the backend is first used, which pins each of PoCL's threads to a core. Tools that pin Julia's threads, like [ThreadPinning.jl](https://github.com/carstenbauer/ThreadPinning.jl), don't affect the threads that run the kernels.
+2) Choose the number of threads with `julia -t N` or the `JULIA_KA_CPU_THREADS` environment variable; see [`CPU`](@ref).
+3) Initialize your data in parallel(!), with a kernel. Under the ["NUMA first-touch policy"](https://queue.acm.org/detail.cfm?id=2513149#:~:text=This%20is%20called%20the%20first,policy%20associated%20with%20a%20task.) a page of memory is placed in the memory domain of the thread that first writes to it. `KernelAbstractions.zeros(backend, dtype, N)` and `fill!` write from the calling thread, which places all the memory in that thread's domain. The example instead allocates with `KernelAbstractions.allocate` and writes the initial values with a kernel (`init = :parallel`), so that the threads that run the computational kernel touch the memory first.
 
 
 **Demonstration:**
 
-If above example is run with 128 Julia threads on a Noctua 2 compute node (128 physical cores distributed over two AMD Milan 7763 CPUs with 4 NUMA domains each), one may get the following numbers (comments for demonstration purposes):
+So far, the example has only been measured with KernelAbstractions 0.10 on a system with a single memory domain, where pinning and the initialization don't change where memory is placed. With 16 Julia threads on an AMD Ryzen 9 5950X (16 physical cores, 1 NUMA domain), one gets the following numbers (comments for demonstration purposes):
 
 ```
-Memory Bandwidth (GB/s): 145.64 # backend = CPU(), init = :parallel
-Compute (GFLOP/s): 24.27
+Memory Bandwidth (GB/s): 41.19 # POCL_AFFINITY=1, init = :parallel
+Compute (GFLOP/s): 6.87
 
-Memory Bandwidth (GB/s): 333.83 # backend = CPU(; static=true), init = :parallel
-Compute (GFLOP/s): 55.64
+Memory Bandwidth (GB/s): 41.22 # POCL_AFFINITY=1, init = :serial
+Compute (GFLOP/s): 6.87
 
-Memory Bandwidth (GB/s): 32.74 # backend = CPU(), init = :serial
-Compute (GFLOP/s): 5.46
+Memory Bandwidth (GB/s): 41.24 # POCL_AFFINITY unset, init = :parallel
+Compute (GFLOP/s): 6.87
 
-Memory Bandwidth (GB/s): 32.46 # backend = CPU(; static=true), init = :serial
-Compute (GFLOP/s): 5.41
+Memory Bandwidth (GB/s): 41.16 # POCL_AFFINITY unset, init = :serial
+Compute (GFLOP/s): 6.86
 ```
 
-The key observations are the following:
-
-* Serial initialization leads to subpar performance (at least a factor of 4.5) independent of the chosen CPU backend. This is a manifestation of remark 3 above.
-* The static CPU backend gives >2x better performance than the one based on the dynamic `Threads.@spawn`. This is a manifestation of remark 2 (and, in some sense, also 1) above.
+As expected for a single memory domain, the four configurations don't differ: the kernel is limited by the memory bandwidth of the system. How much pinning and parallel initialization gain on a system with multiple memory domains has not been measured with the PoCL-based `CPU` backend yet.

docs/src/index.md

@@ -133,6 +133,22 @@ end
 ```
 
 ### 0.10
+- KernelAbstractions requires Julia 1.10 or later.
+- The `CPU` backend is an OpenCL backend: `CPU` is an alias of `POCLBackend`, which compiles
+  kernels like the GPU backends do and runs them with [PoCL](https://portablecl.org) on
+  PoCL's own threads, and no longer on Julia tasks. `CPU(; static=true)` has been removed,
+  as there is no dynamic task scheduling to opt out of anymore.
+- KernelAbstractions is built on [KernelInterface](@ref kernelinterface), which defines
+  `Backend` and the host-side functions (`allocate`, `synchronize`, …) that
+  KernelAbstractions re-exports. User code is unaffected, but backends have to implement
+  KernelInterface; see the [notes for backend implementations](@ref implementations_notes).
+- The Enzyme extension has been removed temporarily, and is planned to return in a later
+  0.10 release.
+- [`KernelAbstractions.@spawn`](@ref) launches kernels from a task like `Threads.@spawn`
+  does, and additionally orders them after the work the spawning task has queued, and can
+  select the device the task uses; see the [Quickstart](@ref).
+- `KernelAbstractions.isgpu` has been removed. Like for `GPU` below, query a capability
+  instead.
 - `ndrange` entries may be index ranges, given statically (`kernel(backend, workgroupsize, (-2:N+3, 0:M+1))`)
   or at launch (`ndrange=(-2:N+3, 0:M+1)`, a single range, or a `CartesianIndices`).
   `@index(Global, Cartesian)` and `@index(Global, NTuple)` return the shifted indices.
@@ -165,6 +181,26 @@ end
   `GPUCompiler.SHARED_METHOD_TABLE`, which backends that override
   `GPUCompiler.method_table_view` have to include, e.g., by implementing
   `GPUCompiler.method_tables` instead.
+- [`@Const`](@ref) arguments are constified on the host when the kernel is launched, instead
+  of at the start of the kernel. Arguments whose Adapt.jl rules don't compile for GPUs, e.g.
+  a `StepRangeLen` or `LinRange`, can now be passed as `@Const`.
+  Backends have to convert kernel arguments with Adapt.jl, both in `KI.argconvert` and in `KI.launch`,
+  since that conversion now constifies them; see the [notes for backend implementations](@ref implementations_notes).
+- Launching a kernel on the `CPU` backend costs more than in 0.9: the kernel is handed to
+  PoCL's threads and the launching task waits for them, which takes about 5 µs per launch
+  even for an empty kernel. Code that launches many small kernels suffers the most. For
+  example, a 7-point stencil over a 4×4×4 range took about 3.5 µs per launch, against 0.2 µs
+  with 0.9. Over a 16×16×16 range, which 0.9 also ran on all threads, both took about 7 µs.
+- On the `CPU` backend, don't pass a fixed workgroup size such as `64`. Omit it
+  (`kernel(CPU())` instead of `kernel(CPU(), 64)`), so that it is automatically chosen for every launch. A
+  one-dimensional workgroup size pads a multidimensional range: with a workgroup size of 64
+  and a range of `(16, 16, 16)`, each of the 256 workgroups has 64 work-items, of which 16
+  are inside the range. A 7-point stencil over that range took about 13 µs per launch with a
+  workgroup size of 64, and 7 µs with the chosen one. Over a 128×128×128 range the chosen
+  size was still 5 to 15% faster.
+- On the `CPU` backend, every work-item checks whether it is inside the range, even when the
+  workgroups cover the range exactly. This can make kernels with a cheap body slightly slower
+  than with 0.9.
 - `foreach_index(f, A)` runs `f` once per index of the array `A` without writing a kernel out,
   and `foreach_index(f, backend, indices)` once per index in a range or `CartesianIndices`.
 

docs/src/quickstart.md

@@ -69,6 +69,36 @@ synchronize(dev)
 
 See also [Memcopy with static NDRange](@ref memcopy_static).
 
+### [Opting out of the bounds check](@id unsafe_indices)
+
+The `ndrange` is padded to a multiple of the workgroup size, so a kernel checks for every
+work-item whether it is inside the `ndrange` before running the body, and the body only
+runs for the work-items that are. A kernel can opt out of the check with `unsafe_indices = true`.
+Its body then also runs for the padded work-items, so it has to test the index itself before touching memory:
+
+```@example mul2_kernel
+@kernel unsafe_indices = true function mul2_unsafe_kernel(A)
+    g = @index(Group, Linear)
+    l = @index(Local, Linear)
+    I = (g - 1) * prod(@groupsize()) + l
+    if I <= length(A)
+        A[I] = 2 * A[I]
+    end
+end
+
+mul2_unsafe_kernel(dev, 64)(A, ndrange = length(A))
+synchronize(dev)
+@assert all(A .== 8)
+```
+
+Don't use `@index(Global)` in any form in such a kernel. Compute the global index from
+`@index(Group)`, `@index(Local)` and [`@groupsize`](@ref) instead, as above, and compare it
+with the size of the range.
+
+!!! note
+    The padded work-items also take part in [`@synchronize`](@ref),
+    so they must reach every barrier that the other work-items of their workgroup reach.
+
 ## Launching kernel on the backend
 
 To launch the kernel on a backend (e.g., `CUDABackend()`, `ROCBackend()`, `oneAPIBackend()`, `MetalBackend()`), we generate the kernel

examples/numa_aware.jl

@@ -1,17 +1,19 @@
 # EXCLUDE FROM TESTING
+# Run with `POCL_AFFINITY=1` to pin the threads that run the kernels, e.g.
+#   POCL_AFFINITY=1 julia -t 128 examples/numa_aware.jl
 using BenchmarkTools
-using Statistics
-using Random
-using ThreadPinning
 using KernelAbstractions
 
-ThreadPinning.pinthreads(:numa)
-
 @kernel function saxpy_kernel(a, @Const(X), Y)
     I = @index(Global)
     @inbounds Y[I] = a * X[I] + Y[I]
 end
 
+@kernel function fill_kernel(A, x)
+    I = @index(Global)
+    @inbounds A[I] = x
+end
+
 """
   measure_membw(; kwargs...) -> membw, flops
 
@@ -24,16 +26,21 @@ function measure_membw(
     )
     bytes = 3 * sizeof(dtype) * N # num bytes transferred in SAXPY
     flops = 2 * N # num flops in SAXY
+    workgroup_size = 1024
 
     a = dtype(3.1415)
+    X = KernelAbstractions.allocate(backend, dtype, N)
+    Y = KernelAbstractions.allocate(backend, dtype, N)
     if init == :serial
-        X = rand(dtype, N)
-        Y = rand(dtype, N)
+        # The calling thread touches all the memory first
+        fill!(X, dtype(1))
+        fill!(Y, dtype(2))
     else
-        X = copyto!(KernelAbstractions.zeros(backend, dtype, N), rand(dtype, N))
-        Y = copyto!(KernelAbstractions.zeros(backend, dtype, N), rand(dtype, N))
+        # The threads that run the kernels touch the memory first
+        fill_kernel(backend, workgroup_size)(X, dtype(1), ndrange = size(X))
+        fill_kernel(backend, workgroup_size)(Y, dtype(2), ndrange = size(Y))
+        KernelAbstractions.synchronize(backend)
     end
-    workgroup_size = 1024
 
     t = @belapsed begin
         kernel = saxpy_kernel($backend, $workgroup_size, $(size(Y)))
@@ -51,10 +58,7 @@ function measure_membw(
     return mem_rate, flop_rate
 end
 
-# Static should be much better (on a system with multiple NUMA domains)
 measure_membw(CPU());
-measure_membw(CPU(; static = true));
 
-# The following has significantly worse performance (even on systems with a single memory domain)!
+# On a system with multiple NUMA domains, this places all the memory in the domain of the calling thread
 # measure_membw(CPU(); init=:serial);
-# measure_membw(CPU(; static=true); init=:serial);