# Cumulative sum on GPUArray using KernelAbstractions

**URL:** <https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098>\
**Category:** GPU\
**Tags:** gpu, gpuarrays, kernelabstractions\
**Created:** [December 22, 2024, 5:12pm UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098 "2024-12-22T17:12:27Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [December 22, 2024, 5:12pm UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098/1 "2024-12-22T17:12:27Z")

</div>

Hello,

I need to implement the cumulative sum (`cumsum`) on a GPU array (CUDA.jl or Metal.jl). Looking at the CUDA.jl repository, I found this definition

```julia
function cumsum!(sums)
    shift = 1

    while shift < length(sums)
        to_add = 0
        @inbounds if threadIdx().x - shift > 0
            to_add = sums[threadIdx().x - shift]
        end

        sync_threads()
        @inbounds if threadIdx().x - shift > 0
            sums[threadIdx().x] += to_add
        end

        sync_threads()
        shift *= 2
    end
end

```

which can be executed inside a CUDA kernel.

Now, I need to implement it using KernelAbstractions.jl, but I’m not familiar with shared memory, especially with KernelAbstractions.jl. The easiest way I thought was something like

```julia
function cumsum!(sums)
    idx = @index(Global)
    shift = 1

    while shift < length(sums)
        to_add = 0
        @inbounds if idx - shift > 0
            to_add = sums[idx - shift]
        end

        KernelAbstractions.@syncronize()
        @inbounds if idx - shift > 0
            sums[idx] += to_add
        end

        KernelAbstractions.@syncronize()
        shift *= 2
    end
end

```

But I don’t know if this is correct, or if I have to use `KernelAbstractions.@localmem` or something else.

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [December 22, 2024, 5:23pm UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098/2 "2024-12-22T17:23:30Z")

</div>

This is a nice [article](https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda) that shows how to implement prefix scan in C/CUDA.

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [December 23, 2024, 9:26am UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098/3 "2024-12-23T09:26:19Z")

</div>

> [@albertomercurio](#):
>
> Now, I need to implement it using [KernelAbstractions.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/KernelAbstractions),

You (probably) don’t have to implement your own kernel for this: CUDA.jl and Metal.jl (as well as AMDGPU.jl) implement `Base.accumulate` for their gpu array types ([CUDA](https://github.com/JuliaGPU/CUDA.jl/blob/master/src/accumulate.jl), [Metal](https://github.com/JuliaGPU/Metal.jl/blob/main/src/accumulate.jl) (and [AMDGPU](https://github.com/JuliaGPU/AMDGPU.jl/blob/master/src/kernels/accumulate.jl) )), so you could just use `accumulate(+, gpu_array)`.

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [December 23, 2024, 9:49am UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098/4 "2024-12-23T09:49:00Z")

</div>

And if you want to port this to KA.jl, the `cumsum!` function you’re looking at is an implementation detail from the sorting kernel, and not something reusable. A more general cumsum kernel can be found in [CUDA.jl/src/accumulate.jl at 972f3f0a9d594c75431df96b027588f8279ad2de · JuliaGPU/CUDA.jl · GitHub](https://github.com/JuliaGPU/CUDA.jl/blob/972f3f0a9d594c75431df96b027588f8279ad2de/src/accumulate.jl), and resembles our mapreduce implementation, which means you can probably draw inspiration from the WIP mapreduce implementation with KA.jl in [Implement mapreduce by vchuravy · Pull Request #561 · JuliaGPU/GPUArrays.jl · GitHub](https://github.com/JuliaGPU/GPUArrays.jl/pull/561).

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [December 24, 2024, 3:16pm UTC](https://discourse.julialang.org/t/cumulative-sum-on-gpuarray-using-kernelabstractions/124098/5 "2024-12-24T15:16:37Z")

</div>

> [@pitsianis](#):
>
> This is a nice [article](https://developer.nvidia.com/gpugems/gpugems3/part-vi-gpu-computing/chapter-39-parallel-prefix-sum-scan-cuda) that shows how to implement prefix scan in C/CUDA.

Thank you, it was very useful.

> [@eldee](#):
>
> You (probably) don’t have to implement your own kernel for this: [CUDA.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/CUDA) and [Metal.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/Metal) (as well as [AMDGPU.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/AMDGPU)) implement `Base.accumulate` for their gpu array types ([CUDA](https://github.com/JuliaGPU/CUDA.jl/blob/master/src/accumulate.jl), [Metal](https://github.com/JuliaGPU/Metal.jl/blob/main/src/accumulate.jl) (and [AMDGPU](https://github.com/JuliaGPU/AMDGPU.jl/blob/master/src/kernels/accumulate.jl) )), so you could just use `accumulate(+, gpu_array)`.

Yes, indeed. I will temporary use their own implementations.

> [@maleadt](#):
>
> A more general cumsum kernel can be found in [CUDA.jl/src/accumulate.jl at 972f3f0a9d594c75431df96b027588f8279ad2de · JuliaGPU/CUDA.jl · GitHub](https://github.com/JuliaGPU/CUDA.jl/blob/972f3f0a9d594c75431df96b027588f8279ad2de/src/accumulate.jl), and resembles our mapreduce implementation, which means you can probably draw inspiration from the WIP mapreduce implementation with [KA.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/KA) in [Implement mapreduce by vchuravy · Pull Request #561 · JuliaGPU/GPUArrays.jl · GitHub](https://github.com/JuliaGPU/GPUArrays.jl/pull/561).

Thanks you, I will give a look.