mapreduce(
    f, op, src, srcs::AbstractArray...;
    backend=nothing,
    init=<none>,
    neutral=nothing,
    acctype=nothing,
    dims=:,
    commutative::Bool=true,
    alg::Algorithm=Auto(),
    workspace=nothing,
)

Apply f to each element of src and reduce the results with the binary operator op. With dims=: (the default; nothing is accepted too) the result is a scalar, returned on the host. With dims an integer or a collection of integers, the result is a new array on src's backend with size 1 along those dimensions; mapreducedim! reduces into an existing array instead.

Both follow the contract that mapreducedim! states: op must be associative, and commutative unless commutative=false, neutral only seeds partial results, and init, when given, is applied once, as op(init, partial). Here the accumulator type, which the result has, is acctype when given, else the type the fold of op settles on from init's type (when given) and the mapped elements: AK.sum(Int8[1, 2]) is an Int, and so is the reduction of a single element, AK.reduce((a, b) -> a + b, [true]) === 1; AK.sum(Int8[1 2]; dims=1, init=Int16(0)) is an array of Ints. An empty reduction gives init as it is, and without init is an ArgumentError: of a whole array, and along dims of any output whose slice is empty (sum, prod and count give zero or one instead). These rules are AcceleratedKernels' own, and differ from Base's in a few places; see Differences from Base.

Several source arrays with equal axes are reduced without materializing the mapped array (f takes one argument per array); mismatched axes throw a DimensionMismatch. A Broadcasted object is also accepted as the single source. Before Julia 1.12, both are materialized when the reduction runs. backend is derived from all sources.

alg is Auto() by default: CPUThreads.Partitioned on the host and BlockReduce on GPUs, with the device's settings. On the host, multithreading reductions only improves performance for operations that hide the memory latency and thread launch overhead, e.g. reductions of tuples or structs, or expensive operators.

workspace takes the scratch memory of a workspace made for the same call, so that the reduction allocates none of its own (it still allocates its result along dims, and before Julia 1.12 a materialized source).

Examples

Computing a sum of squares, reducing down to a scalar that is copied to host:

import AcceleratedKernels as AKusing CUDAv = CuArray{Int16}(rand(1:1000, 100_000))vsumsq = AK.mapreduce(x -> Int(x) * x, +, v)

Computing the sums of squares of the columns of a matrix:

import AcceleratedKernels as AKusing Metalf(x) = x * xm = MtlArray(rand(Int32(1):Int32(100), 10, 100_000))colsumsq = AK.mapreduce(f, +, m; dims=1)        # a 1×100_000 MtlArray

Computing a two-input reduction along dims:

colsums = AK.mapreduce((x, y) -> x * y, +, a, b; dims=1)

mapreducedim!(
    f, op, R::AbstractArray, src;
    backend=nothing,
    init=<none>,
    neutral=nothing,
    overwrite::Bool=false,
    acctype=nothing,
    commutative::Bool=true,
    alg::Algorithm=Auto(),
    workspace=nothing,
) -> R

Reduce src (an array or a Broadcasted object) into R: the dimensions along which R has size 1 and src does not are reduced, and the others must match (R may leave off trailing dimensions of size 1). R must not alias src. backend is derived from R and src.

This is the contract of every AcceleratedKernels reduction, modelled on CUB's rather than on Base's (see Differences from Base).

op must be associative, and by default commutative: as CUB's reductions do, AK then combines elements in an order that depends on the algorithm, its settings and the shape, not in element order. With commutative=false, elements are combined in their order (each partial result reduces consecutive elements, and only neighbouring partial results are combined, the earlier on the left), so that operators such as matrix products are reduced correctly; the grouping is still unspecified, as in Base. The GPU kernels that keep the order can be slower, mostly for inputs that fit in the device's cache. Either way the order is fixed for a given setting and shape, so results are reproducible on one device.

Partial results start from neutral, which must be a two-sided identity of op and never appears in the result. It defaults to GPUArraysCore.neutral_element(op, T) where that is defined; otherwise each partial result starts from its first element, Base.mapreduce_first(f, op, x), so no neutral element is needed.

Partial results have one type, the accumulator type, and are converted to eltype(R) only when stored. It is acctype when that is given; otherwise the type the fold of op settles on, starting from eltype(R) and the mapped elements. Every element is also a one-element partial result, and partial results are combined with each other, so those types are joined in as well (R's values and init are applied once, and are not partial results). A sum of Int8s accumulates in Int, and Float32s reduced into a Float64 array accumulate in Float64. An acctype that cannot hold the partial results at all is an ArgumentError; whether their values fit is the caller's obligation (a value that does not fit throws an InexactError, where the backend reports errors thrown in kernels).

For each output whose slice is not empty, with partial its reduction, R[i] becomes op(init, partial) when init is given (applied once); else partial with overwrite=true; else op(R[i], partial), folding into R's previous value as Base.mapreducedim! does. An output whose slice is empty is set to init when that is given, and is otherwise not written, with or without overwrite.

Mismatched shapes are a DimensionMismatch. An R that aliases src, an op that inference shows always throws for these types (when there is something to combine; an explicit acctype that op cannot combine at all is rejected whatever the input), a non-bits accumulator type for a kernel algorithm, and an algorithm that cannot run on the backend are ArgumentErrors.

workspace takes the scratch memory of a workspace made for the same call.

import AcceleratedKernels as AKusing CUDAA = CuArray([1 3; 2 4])R = CuArray([10 20])AK.mapreducedim!(identity, +, R, A)                  # [13 27]AK.mapreducedim!(identity, +, R, A; overwrite=true)  # [3 7]AK.mapreducedim!(identity, +, R, A; init=100)        # [103 107]