maleadt · Sep 30, 2026 · Sep 23, 2026 · Sep 23, 2026 · Sep 23, 2026 · Sep 24, 2026
diff --git a/.github/workflows/CI-CPU.yml b/.github/workflows/CI-CPU.yml
@@ -102,26 +102,37 @@ jobs:
       - name: Promote OpenCL to test dependency
         run: julia test/promote.jl opencl
       - uses: julia-actions/julia-runtest@v1
-  # cpuKA:
-  #   name: KA CPU Backend
-  #   runs-on: ubuntu-latest
-  #   timeout-minutes: 60
-  #   permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
-  #     actions: write
-  #     contents: read
-  #   strategy:
-  #     fail-fast: true
-  #   steps:
-  #     - uses: actions/checkout@v7
-  #     - uses: julia-actions/setup-julia@v3
-  #       with:
-  #         version: 1
-  #         arch: x64
-  #     - uses: julia-actions/cache@v3
-  #     - uses: julia-actions/julia-buildpkg@v1
-  #     - uses: julia-actions/julia-runtest@v1
-  #       with:
-  #         test_args: '--cpu-ka'
+  cpuKA:
+    # KernelAbstractions 0.10 (unreleased, `main`) runs AcceleratedKernels' kernels on the host
+    # backend, through PoCL; `--cpu-ka` tests them there.
+    name: KA 0.10 host kernels - Julia ${{ matrix.version }}
+    runs-on: ubuntu-latest
+    timeout-minutes: 60
+    permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created
+      actions: write
+      contents: read
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.10'
+          - '1.13'
+    steps:
+      - uses: actions/checkout@v7
+      - uses: julia-actions/setup-julia@v3
+        with:
+          version: ${{ matrix.version }}
+          arch: x64
+      - uses: julia-actions/cache@v3
+      - name: Develop AcceleratedKernels and KernelAbstractions main into the test environment
+        run: |
+          julia --project=test -e '
+            using Pkg
+            Pkg.develop([PackageSpec(path="."),
+                         PackageSpec(url="https://github.com/JuliaGPU/KernelAbstractions.jl")])'
+      - name: Run tests
+        # only the generic tests: the others, Aqua included, do not depend on the backend
+        run: julia --project=test test/runtests.jl --cpu-ka cpu-ka/
   docs:
     name: Documentation
     runs-on: ubuntu-latest

diff --git a/Project.toml b/Project.toml
@@ -32,9 +32,9 @@ AMDGPU = "1.3.4, 2"
 ArgCheck = "2"
 Atomix = "0.1, 1"
 CUDA = "6"
-CUDACore = "6"
-GPUArraysCore = "0.2.0"
-KernelAbstractions = "0.9.34, 0.10"
+CUDACore = "6.4.1"
+GPUArraysCore = "0.2.1"
+KernelAbstractions = "0.9.43, 0.10"
 Markdown = "1"
 Metal = "1.10.1"
 OpenCL = "0.10"

diff --git a/README.md b/README.md
@@ -270,12 +270,8 @@ If you need other algorithms in your work that may be of general use, please ope
 | [General Looping](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/foreachindex/) | `foreachindex`, `foraxes`                        | `Kokkos::parallel_for` `RAJA::forall` `thrust::transform` |
 | [Mapping](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/map/) | `map` `map!`                                     | `thrust::transform`                                       |
 | [Sorting](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/sort/) | `sort` `sort!`                                   | `sort` `sort_team` `stable_sort`                          |
-|                                               | `sample_sort!` `sample_sortperm!`                |                                                           |
-|                                               | `merge_sort` `merge_sort!`                       |                                                           |
-|                                               | `merge_sort_by_key` `merge_sort_by_key!`         | `sort_team_by_key`                                        |
 |                                               | `sortperm` `sortperm!`                           | `sort_permutation` `index_permutation`                    |
-|                                               | `merge_sortperm` `merge_sortperm!`               |                                                           |
-|                                               | `merge_sortperm_lowmem` `merge_sortperm_lowmem!` |                                                           |
+|                                               | `sort_by_key!`                                   | `sort_by_key` `sort_team_by_key` `cub::DeviceRadixSort::SortPairs` |
 | [Reduction](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/reduce/) | `reduce`                                         | `Kokkos:parallel_reduce` `fold` `aggregate`               |
 | [MapReduce](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/mapreduce/) | `mapreduce`                                      | `transform_reduce` `fold`                                 |
 | [Accumulation](https://juliagpu.github.io/AcceleratedKernels.jl/stable/api/accumulate/) | `accumulate` `accumulate!`                       | `prefix_sum` `thrust::scan` `cumsum`                      |
@@ -345,6 +341,12 @@ Start Julia with multiple threads to run the tests on a multithreaded CPU backen
 $> julia --threads=4 -e 'import Pkg; Pkg.test("AcceleratedKernels.jl")'
 ```
 
+With KernelAbstractions 0.10, whose host backend runs kernels on PoCL, `--cpu-ka` tests
+AcceleratedKernels' GPU kernels on host arrays instead of the threaded CPU algorithms:
+```bash
+$> julia -e 'import Pkg; Pkg.test("AcceleratedKernels"; test_args=["--cpu-ka"])'
+```
+
 
 ## 8. Issues and Debugging
 As the compilation pipeline of GPU kernels is different to that of base Julia, error messages also look different - for example, where Julia would insert an exception when a variable name was not defined (e.g. we had a typo), a GPU kernel throwing exceptions cannot be compiled and instead you'll see some cascading errors like `"[...] compiling [...] resulted in invalid LLVM IR"` caused by `"Reason: unsupported use of an undefined name"` resulting in `"Reason: unsupported dynamic function invocation"`, etc.

diff --git a/docs/make.jl b/docs/make.jl
@@ -21,6 +21,7 @@ makedocs(;
         "Performance Tips" => "performance.md",
         "Manual" =>[
             "Using Different Backends" => "api/using_backends.md",
+            "Algorithms and Backends" => "api/algorithms.md",
             "General Loops" => "api/foreachindex.md",
             "Map" => "api/map.md",
             "Sorting" => "api/sort.md",
@@ -32,11 +33,14 @@ makedocs(;
             "Reverse" => "api/reverse.md",
             "Predicates" => "api/predicates.md",
             "Arithmetics" => "api/arithmetics.md",
+            "Scratch Memory" => "api/workspace.md",
             "Custom Structs" => "api/custom_structs.md",
             "Task Partitioning" => "api/task_partition.md",
             "Utilities" => "api/utilities.md",
+            "Differences from Base" => "api/differences.md",
         ],
         "Testing" => "testing.md",
+        "Tuning and Capabilities" => "tuning.md",
         "Debugging Kernels" => "debugging.md",
         "Roadmap" => "roadmap.md",
         "References" => "references.md",

diff --git a/docs/src/api/accumulate.md b/docs/src/api/accumulate.md
@@ -5,3 +5,24 @@
 AcceleratedKernels.accumulate!
 AcceleratedKernels.accumulate
 ```
+
+[`Auto()`](@ref AcceleratedKernels.Auto), the default `alg`, uses
+[`CPUThreads.Partitioned`](@ref AcceleratedKernels.CPUThreads.Partitioned) on the host; on GPUs,
+[`ScanPrefixes`](@ref AcceleratedKernels.ScanPrefixes) for whole arrays and
+[`SliceScan`](@ref AcceleratedKernels.SliceScan) along `dims`, with the device's settings. Pass an
+algorithm to choose it and its settings yourself:
+
+```@docs
+AcceleratedKernels.ScanAlgorithm
+AcceleratedKernels.ScanPrefixes
+AcceleratedKernels.DecoupledLookback
+AcceleratedKernels.SliceScan
+```
+
+```julia
+v = CuArray(rand(Int32(1):Int32(100), 1_000_000))
+AK.accumulate!(+, v; init=Int32(0), alg=AK.ScanPrefixes(block_size=512))  # 512 threads, 8 items per thread
+AK.accumulate!(+, v; init=Int32(0), alg=AK.DecoupledLookback())            # CUDA and AMDGPU only
+m = CuArray(rand(Int32(1):Int32(100), 100, 10_000))
+AK.accumulate(+, m; init=Int32(0), dims=2, alg=AK.SliceScan(block_size=128))
+```
diff --git a/docs/src/api/algorithms.md b/docs/src/api/algorithms.md
@@ -0,0 +1,57 @@
+### Algorithms and backends
+
+Every algorithmic operation takes an `alg` keyword. Its default,
+[`Auto()`](@ref AcceleratedKernels.Auto), lets AcceleratedKernels choose the algorithm and its
+settings; any other value is an algorithm you choose, with the settings you give it.
+
+```@docs
+AcceleratedKernels.Algorithm
+AcceleratedKernels.Auto
+AcceleratedKernels.SortAlgorithm
+AcceleratedKernels.CPUThreads
+```
+
+The algorithm types of each operation are documented with it, for [sorting](sort.md),
+[reductions](reduce.md), [scans](accumulate.md), [`findall`](findall.md) and
+[`any`/`all`](predicates.md).
+
+#### How `Auto` chooses
+
+`Auto` looks at the backend and its current device, the element type, the layout (a whole array or
+slices along `dims`, and their length) and the ordering, never at the array's contents. On the
+host backend it chooses the family's `CPUThreads` algorithm, which runs on Julia threads. On a GPU
+it follows per-device tuning values: for sorting, `BitonicSort` for short arrays and slices,
+`RadixSort` for long whole arrays of the element types it supports, and `MergeSort` otherwise,
+subject to `Auto`'s requirements (`stable=true` by default); for reductions, `BlockReduce`; for
+scans, `ScanPrefixes` for whole arrays and `SliceScan` along `dims`; for `findall`,
+`ScanScatter`; for `any`/`all`, `ConcurrentWrite` (`ViaReduce` on oneAPI).
+
+Algorithms carry their settings as fields, and a field left at `nothing` takes the device's
+tuned value:
+```julia
+AK.sort!(v; alg=AK.RadixSort())                   # the device's block size and items per thread
+AK.sort!(v; alg=AK.RadixSort(block_size=512))     # 512 threads, the device's items per thread
+```
+
+Explicit algorithms are checked before any data is touched: an algorithm the operation, the
+element type, the ordering or the backend cannot run, or an invalid setting, is an
+`ArgumentError`. Each new combination of settings compiles new kernels.
+
+#### The `backend` keyword
+
+Operations run on the backend of their arrays: `backend` is derived from every array argument,
+the destination first, and all of them must agree. Ranges, `CartesianIndices`, `LinearIndices`
+(also through views and reshapes), numbers and other non-array arguments do not count. If no
+argument determines the backend (e.g. a loop over a range), the operation runs on the host.
+
+Pass `backend` explicitly for arrays that cannot tell (a range to be processed on a GPU, for
+instance), or for memory that several backends can access. AcceleratedKernels does not check that
+the backend can reach the arrays. The device and stream are those of the calling task, as set by
+the backend package (e.g. `CUDA.device!`); make sure the arrays live on that device.
+
+#### The host backend
+
+Arrays in host memory (`Array` and views of it) are on the host backend. `Auto` processes them on
+Julia threads with the `CPUThreads` algorithms. On KernelAbstractions 0.10 the host backend also
+runs AcceleratedKernels' GPU kernels (on PoCL), so an explicit kernel algorithm such as
+`MergeSort()` works on host arrays there; on KernelAbstractions 0.9 it is an `ArgumentError`.
diff --git a/docs/src/api/binarysearch.md b/docs/src/api/binarysearch.md
@@ -1,9 +1,9 @@
 ### Binary Search
 
-Find the indices where some elements `x` should be inserted into a sorted sequence `v` to maintain the sorted order. Effectively applying the Julia.Base functions in parallel on a GPU using `foreachindex`.
-- `searchsortedfirst!` (in-place), `searchsortedfirst` (allocating): index of first element in `v` >= `x[j]`.
-- `searchsortedlast!`, `searchsortedlast`: index of last element in `v` <= `x[j]`.
-- **Other names**: `thrust::upper_bound`, `std::lower_bound`.
+Find the indices where many elements `x[j]` would be inserted into a sorted sequence `v` to maintain the sorted order: the Julia Base functions applied to each query in parallel. There are no allocating forms, because `Base.searchsortedfirst(v, x)` treats a vector `x` as one value.
+- `searchsortedfirst!`: index of the first element of `v` not ordered before `x[j]` (`>= x[j]` by default).
+- `searchsortedlast!`: index of the last element of `v` not ordered after `x[j]` (`<= x[j]` by default).
+- **Other names**: `thrust::lower_bound`, `thrust::upper_bound`, `std::lower_bound`.
 
 
 Example:
@@ -13,7 +13,7 @@ using Metal
 
 # Sorted array
 v = MtlArray(rand(Float32, 100_000))
-AK.sort!(v; alg=AK.MergeSort())
+AK.sort!(v)
 
 # Elements `x` to place within `v` at indices `ix`
 x = MtlArray(rand(Float32, 10_000))
@@ -25,7 +25,5 @@ AK.searchsortedfirst!(ix, v, x)
 
 ```@docs
 AcceleratedKernels.searchsortedfirst!
-AcceleratedKernels.searchsortedfirst
 AcceleratedKernels.searchsortedlast!
-AcceleratedKernels.searchsortedlast
 ```
diff --git a/docs/src/api/custom_structs.md b/docs/src/api/custom_structs.md
@@ -14,7 +14,7 @@ using CUDA
 
 function complex_any(x, y)
     # Calling `any` on a normal Julia range, but running on x's backend
-    AK.any(1:length(x), AK.get_backend(x)) do i
+    AK.any(1:length(x); backend=AK.get_backend(x)) do i
         x[i] < 0 && y[i] > 0
     end
 end

diff --git a/docs/src/api/differences.md b/docs/src/api/differences.md
@@ -0,0 +1,25 @@
+### Differences from Base
+
+AcceleratedKernels' functions with Base's names take the arguments Base's do, but follow
+contracts of their own, modelled on GPU libraries such as CUB: `init` is applied once, empty
+inputs have no implicit result, and results have one documented type. Base's rules that exist for
+Base's own reasons (empty results chosen per operator, `typeof(init)` as the result type along
+`dims`, and so on) are left to front-ends: GPUArrays.jl, for example, implements Base's API for
+GPU arrays on top of AcceleratedKernels and reproduces Base's results. The differences:
+
+| Call | Base | AcceleratedKernels |
+|---|---|---|
+| `reduce(op, A)`, `mapreduce(f, op, A)` of an empty `A`, no `init` | `Base.mapreduce_empty`: e.g. `sum(Int[]) == 0`, an error for `maximum` and for most mapped reductions | `ArgumentError` (`sum` and `prod` give zero and one) |
+| `mapreduce(f, op, A; dims)` with an empty reduced dimension, no `init` | `Base.reducedim_init`: e.g. `[0 0]` for `x -> x + 1` with `+`, an error for `max` | `ArgumentError` where there are outputs (`sum`, `prod` and `count` give zero or one) |
+| a one-element reduction, e.g. `reduce((a, b) -> a + b, [true])` | `true` (`Base.mapreduce_first`) | `1`, the accumulator type |
+| the element type of `mapreduce(f, op, A; dims, init)` | `typeof(init)` | the accumulator type: `sum(Int8[1 2]; dims=1, init=Int16(0))` is a `Matrix{Int}` |
+| the accumulator type of reductions into an array (`sum!`, reductions along `dims`) | depends on the code path, e.g. on the reduced dimension | one rule, the fold type from `eltype(R)` (see [`mapreducedim!`](@ref AcceleratedKernels.mapreducedim!)) |
+| an `init` that is not a neutral element of `op` | outside the contract: `init` must be neutral, and it is unspecified whether it is used for non-empty collections | any value, applied exactly once |
+| `reduce(op, A; dims)` for an `op` without `Base.reducedim_init`, such as a closure | `MethodError` | works |
+| a non-commutative `op` in a reduction | works (elements keep their order) | unsupported |
+| the running-value type of `accumulate!(op, B, A)` | depends on the code path: the fold type from `init` or the elements for vectors and `dims=1`, `eltype(B)` along `dims ≥ 2` | the fold type from `eltype(B)` joined with `init`'s type and the elements (unless `acctype` is given), so for the usual operators not narrower than `eltype(B)` |
+| the element type of `accumulate(op, A)` on Julia 1.10 | `Base.promote_op(op, T, T)` | the fold type, as Base's on Julia 1.13 |
+| `accumulate(op, A; dims, init)` with `dims > ndims(A)` | copies `A`, ignoring `init` | applies `init` to every element |
+| a non-associative `op` in a scan, such as `-` | works (a sequential recurrence) | unsupported |
+| `findall(pred, A)` of a 0-dimensional `A` | `Int` indices | `CartesianIndex{0}`, `keys(A)`'s, unless `items=LinearIndices(A)` |
+| `any`, `all` with a predicate that returns `missing` | three-valued logic | an error: the predicate must return a `Bool` (an `ArgumentError` before launching where inference shows it cannot) |
diff --git a/docs/src/api/findall.md b/docs/src/api/findall.md
@@ -2,5 +2,22 @@
 
 ```@docs
 AcceleratedKernels.findall
+```
+
+`findall` selects items: indices by default, or anything else of the input's length, such as
+the values themselves.
+
+```julia
+v = CuArray(rand(Float32, 1000))
+AK.findall(x -> x > 0.5f0, v)                  # indices
+AK.findall(x -> x > 0.5f0, v; items=v)         # the selected values, as v[v .> 0.5f0]
+```
+
+[`Auto()`](@ref AcceleratedKernels.Auto) uses
+[`CPUThreads.Partitioned`](@ref AcceleratedKernels.CPUThreads.Partitioned) on the host and
+[`ScanScatter`](@ref AcceleratedKernels.ScanScatter) on GPUs, with the device's settings.
+
+```@docs
+AcceleratedKernels.FindallAlgorithm
 AcceleratedKernels.ScanScatter
 ```
diff --git a/docs/src/api/foreachindex.md b/docs/src/api/foreachindex.md
@@ -56,14 +56,22 @@ v2 = oneArray(rand(Float32, 100_000))
 f(v1, v2)
 ```
 
+The loop runs on the backend of the iterable. A range does not determine one, so pass `backend` when a loop over a range accesses GPU arrays:
+```julia
+AK.foreachindex(1:n; backend=AK.get_backend(v)) do i
+    v[i] = i
+end
+```
+Host arrays and ranges run on Julia threads.
+
 All GPU functions allow you to specify a block size - this is often a power of two (mostly 64, 128, 256, 512); the optimum depends on the algorithm, input data and hardware - you can try the different values and `@time` or `@benchmark` them:
 ```julia
-@time AK.foreachindex(f, itr_gpu, block_size=512)
+@time AK.foreachindex(f, itr_gpu; block_size=512)
 ```
 
 Similarly, for performance on the CPU the overhead of spawning threads should be masked by processing more elements per thread (but there is no reason here to launch more threads than `Threads.nthreads()`, the number of threads Julia was started with); the optimum depends on how expensive `f` is - again, benchmarking is your friend:
 ```julia
-@time AK.foreachindex(f, itr_cpu, max_tasks=16, min_elems=1000)
+@time AK.foreachindex(f, itr_cpu; max_tasks=16, min_elems=1000)
 ```
 
 

diff --git a/docs/src/api/mapreduce.md b/docs/src/api/mapreduce.md
@@ -8,3 +8,12 @@ Equivalent to `reduce(op, map(f, iterable))`, without saving the intermediate ma
 ```@docs
 AcceleratedKernels.mapreduce
 ```
+
+To reduce into an existing array, for example to accumulate into it across calls or to avoid an
+allocation, use `mapreducedim!`. Its docstring states the contract every reduction follows:
+the operator algebra, neutral elements, the accumulator type, how `init` and the destination's
+values combine with the result, and empty slices.
+
+```@docs
+AcceleratedKernels.mapreducedim!
+```
diff --git a/docs/src/api/predicates.md b/docs/src/api/predicates.md
@@ -9,4 +9,18 @@ AcceleratedKernels.any
 AcceleratedKernels.all
 ```
 
-**Note on the `cooperative` keyword**: some older platforms crash when multiple threads write to the same memory location in a global array (e.g. old Intel Graphics); if all threads were to write the same value, it is well-defined on others (e.g. CUDA F4.2 says "If a non-atomic instruction executed by a warp writes to the same location in global memory for more than one of the threads of the warp, only one thread performs a write and which thread does it is undefined."). This "cooperative" thread behaviour allows for a faster implementation; if you have a platform - the only one I know is Intel UHD Graphics - that crashes, set `cooperative=false` to use a safer `mapreduce`-based implementation.
+[`Auto()`](@ref AcceleratedKernels.Auto) uses
+[`CPUThreads.Partitioned`](@ref AcceleratedKernels.CPUThreads.Partitioned) on the host and
+[`ConcurrentWrite`](@ref AcceleratedKernels.ConcurrentWrite) on GPUs, in which many threads write
+the same value to one memory location. That is well-defined (CUDA F4.2: "If a non-atomic
+instruction executed by a warp writes to the same location in global memory for more than one of
+the threads of the warp, only one thread performs a write and which thread does it is
+undefined."), but some older platforms (Intel UHD Graphics) have been reported to hang on it, so
+on oneAPI `Auto` uses the `mapreduce`-based [`ViaReduce`](@ref AcceleratedKernels.ViaReduce)
+instead. An explicit `alg=ConcurrentWrite()` runs on every GPU.
+
+```@docs
+AcceleratedKernels.PredicateAlgorithm
+AcceleratedKernels.ConcurrentWrite
+AcceleratedKernels.ViaReduce
+```