Accumulate / Prefix Sum / Scan
AcceleratedKernels.accumulate! — Function
accumulate!(op, v::AbstractArray; kwargs...) -> v
accumulate!(op, dst::AbstractArray, src::AbstractArray; kwargs...) -> dst
accumulate!(
op, dst, src;
backend=nothing,
init=<none>,
neutral=nothing,
dims::Union{Nothing, Integer}=nothing,
inclusive::Bool=true,
acctype=nothing,
alg::Algorithm=Auto(),
workspace=nothing,
)Compute accumulated running totals along a sequence by applying a binary operator to all elements up to the current one; often used in GPU programming as a first step in finding / extracting subsets of data. The first form scans v in place; the second writes the scan of src to dst, which must be src itself or not overlap it.
Also known as a prefix sum, thrust::scan or cub::DeviceScan; it is inclusive (or exclusive) if the first element is included in the accumulation (or not).
With dims=nothing the array is scanned in linear order, whatever its number of dimensions; with an integer dims, each slice along that dimension is scanned independently. A dims beyond the array's dimensions makes every slice one element long.
This contract is AcceleratedKernels' own (see Differences from Base).
op must be associative, as elements are combined in parallel; it need not be commutative, as every combination keeps the elements in order (matrix products are fine). Non-associative operators such as - are not supported.
An inclusive scan (the default) computes the running values y[1] = Base.reduce_first(op, src[1]), or op(init, src[1]) with init, then y[k] = op(y[k-1], src[k]), and stores dst[k] = y[k]; init is applied once per slice. An exclusive scan (inclusive=false) starts from y[1] = init, or from the neutral element of op without it (an ArgumentError when none is known), and continues with y[k] = op(y[k-1], src[k-1]).
Partial results start from neutral, a two-sided identity of op. It defaults to GPUArraysCore.neutral_element(op, T) where that is defined; for other operators, partial results start from their first element instead, so inclusive scans need no neutral element.
The running values have one type, the accumulator type, and are converted to eltype(dst) only when stored. It is acctype when that is given; otherwise the type the fold of op settles on, starting from eltype(dst) (joined with init's type) and the elements, as for mapreducedim!: Int for add_sum over Int8, and Float64 for Float32s scanned into a Float64 array. Where it differs from eltype(dst), the scan runs in a scratch array of that type. An acctype that cannot hold the running values at all is an ArgumentError; whether their values fit is the caller's obligation. accumulate allocates the fold type from init's type and the elements.
alg is Auto() by default: CPUThreads.Partitioned on the host, and on GPUs ScanPrefixes for whole arrays and SliceScan along dims, with the device's settings. backend is derived from dst and src.
On the host, accumulation is typically a memory-bound operation, so multithreaded accumulation only becomes faster for more compute-heavy operations that hide memory latency, e.g. accumulating tuples or structs, or expensive operators.
workspace takes the scratch memory of a workspace made for the same call, so that the scan allocates none of its own.
Examples
Example computing an inclusive prefix sum (the typical GPU "scan"):
import AcceleratedKernels as AKusing oneAPIv = oneAPI.ones(Int32, 100_000)AK.accumulate!(+, v) # 1, 2, 3, ...AK.accumulate!(+, v; inclusive=false) # 0, 1, 3, ...: starts from the neutral element# Choose the algorithm and its settingsAK.accumulate!(+, v; alg=AK.ScanPrefixes(block_size=512))AcceleratedKernels.accumulate — Function
accumulate(op, v::AbstractArray; init=<none>, kwargs...)Out-of-place version of accumulate!, with the same keywords. The result's element type is acctype when it is given, else the type the fold of op settles on from init's type (when given) and the elements, e.g. an Int array for accumulate(+, Int8[1, 2]; init=0) and for accumulate(Base.add_sum, Int8[1, 2]).
Auto(), the default alg, uses CPUThreads.Partitioned on the host; on GPUs, ScanPrefixes for whole arrays and SliceScan along dims, with the device's settings. Pass an algorithm to choose it and its settings yourself:
AcceleratedKernels.ScanAlgorithm — Type
ScanAlgorithm <: AlgorithmSupertype of the scan (accumulate) algorithms: ScanPrefixes and SliceScan. Scans also accept CPUThreads.Partitioned.
AcceleratedKernels.ScanPrefixes — Type
ScanPrefixes(; block_size=nothing, items_per_thread=nothing)GPU scan of a whole array in tiles of block_size * items_per_thread elements (block_size a power of two up to 1024): each block reduces its tile, the tiles' totals are scanned the same way, and each block then scans its tile from the total of the tiles before it. On devices whose ScanTuning says so (CPUs), the blocks instead scan their tile first and then add the total of the tiles before it. No block waits for another, so it runs on every backend. The default items_per_thread is the largest, up to 8, whose tile fits in the tuning's local-memory budget.
AcceleratedKernels.SliceScan — Type
SliceScan(; block_size=nothing, items_per_thread=nothing)GPU scan of each slice along dims. Short slices, and slices whose neighbours are next to them in memory (when the dimensions before dims hold several elements), are scanned by one thread each; other slices by one block each, in tiles of block_size * items_per_thread consecutive elements (block_size a power of two up to 1024, smaller for short slices). When that leaves too few threads to fill the device, each slice is split into chunks, which takes two more passes: one over the source to reduce each chunk, and a scan of the chunks' totals.
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 threadm = CuArray(rand(Int32(1):Int32(100), 100, 10_000))AK.accumulate(+, m; init=Int32(0), dims=2, alg=AK.SliceScan(block_size=128))