Arrays, on CPU and GPU

Array-of-structs is the default

Array{Double32} stores each element as its 8-byte pair, inline. It is the right representation for elementwise work on both CPU and GPU: one load per element, ordinary broadcast fusion, ordinary similar/copyto!/views.

julia> v = Double32[1, 2, 3];

julia> v .+= one(Double32) ./ 3;

julia> Float64.(v)
3-element Vector{Float64}:
 1.3333333333333321
 2.3333333333333357
 3.3333333333333357

On the GPU this is CuArray{Double32} / ROCArray{Double32} — Double32 is isbits, so device arrays hold it inline and kernels receive it in registers:

using CUDA
A = CuArray(rand(Double32, 10^6))
B = A .* A .- Double32(0.5f0)              # fused device kernel

Reductions through the Base generics

The important patterns require no package function. Once the scalar methods exist, Base's mapreduce machinery supplies the rest, on the CPU and through GPUArrays on a device:

julia> A32 = Float32[1f8; fill(1.0f0, 1000)];   # big value then many small ones

julia> Float32(sum(A32))                        # Float32 accumulator: wrong tail
1.0000097f8

julia> Float64(sum(Double32, A32))              # Double32 accumulator: exact
1.00001e8

Package names exist only where Base has no spelling:

julia> x32 = Float32[1, 2, 3]; y32 = Float32[4, 5, 6];

julia> Double32s.compensated_dot(x32, y32)      # each product formed EXACTLY
32.0

Double32s.compensated_dot forms every Float32×Float32 product exactly as a pair (two_prod_) before accumulating, so no information is lost until the accumulation itself — strictly better than rounding each product to Float32 first. Double32s.compensated_norm2 adds an exact power-of-two scaling so squares neither overflow nor underflow:

julia> Double32s.compensated_norm2(Double32[1f30, 1f30])   # Float32 squares would overflow
1.414213583653427e30

Double32s.compensated_sum and Double32s.compensated_mapreduce exist for one reason: they pin a sequential order for reproducibility work. Prefer sum(Double32, A) otherwise — it is the parallel AcceleratedKernels tree, reproducible per backend for a pinned block size.

Struct-of-arrays: SplitDouble32Array

Some device kernels want the two words as independent coalesced streams:

• a seed phase reads only high words (Float32 speed, half the bandwidth);
• a vendor Float32 routine produces the high approximation, a correction pass fills in the low words;
• a reduction keeps high and low partials in separate scratch arrays.

julia> v = [Double32(1.0f0) / 3, Double32(2.0f0)];

julia> S = Double32s.split(v);                  # two Float32 arrays

julia> S.hi, S.lo
(Float32[0.33333334, 2.0], Float32[-9.934108f-9, 0.0])

julia> S[1] === v[1]                            # still an AbstractArray{Double32}
true

julia> S[2] = Double32(5.0f0, 2.0f0^-28); S.hi[2]
5.0f0

julia> Double32s.unsplit(S) == [v[1], Double32(5.0f0, 2.0f0^-28)]
true

SplitDouble32Array is a full AbstractArray{Double32}: indexing reconstructs the pair, setindex! stores both words, and both backing arrays must share axes — the constructor enforces it. Its indexing methods propagate @inbounds, which matters inside GPU kernels where a bounds check drags in exception machinery.

Constructing a split array on a device

There are three routes, all of which keep the data resident on the device.

Split an existing array. split and unsplit are broadcasts, so they stay on whatever backend owns the data, with no host round trip and no scalar indexing:

using CUDA
A = CuArray(rand(Double32, 10^6))

S = Double32s.split(A)                      # two CuArray{Float32}
S.hi                                        # Float32-only pass reads this alone
A2 = Double32s.unsplit(S)                   # back to CuArray{Double32}

Allocate for a backend. Double32s.allocate_split needs KernelAbstractions loaded (CUDA and AMDGPU both bring it):

using KernelAbstractions, CUDA
S = Double32s.allocate_split(CUDABackend(), 10^6)
KernelAbstractions.get_backend(S)           # CUDABackend()

Wrap arrays you built yourself, for interop with a vendor routine that produced the two streams separately:

S = Double32s.SplitDouble32Array(CuArray(h_hi), CuArray(h_lo))
seed = some_float32_vendor_pass(S.hi)       # Float32-only first pass
# ... correction kernel then consumes S.hi, S.lo, writes a Double32 result

Either way it goes into kernels directly: Adapt recurses into both words, so the wrapper is rebuilt around device arrays automatically.

@kernel function scale!(S, f)
    i = @index(Global)
    @inbounds S[i] = S[i] * f
end

scale!(CUDABackend())(S, Double32(2.0f0); ndrange = length(S))

Choosing a representation

The design's measured recommendation (section 13) is array-of-structs as the primary user representation, and structure-of-arrays where a kernel genuinely consumes the two words separately. The table is a starting point; where a measurement on your own workload disagrees with it, the measurement governs.