# Correct implementation of CuArray's slicing operations

**URL:** <https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600>\
**Category:** GPU\
**Created:** [November 21, 2022, 6:39pm UTC](https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600 "2022-11-21T18:39:46Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![FujiwaraTakumiEH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fujiwaratakumieh/32/37975_2.png) [@FujiwaraTakumiEH](https://discourse.julialang.org/u/FujiwaraTakumiEH)\
**Post date:** [November 21, 2022, 6:39pm UTC](https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600/1 "2022-11-21T18:39:46Z")

</div>

Hello, sometimes I need slicing operations on CuArray. I would like to know the correct (and efficient) way to achieve this. I made 4 benchmarks:

## Datasets:

```julia
using BenchmarkTools
using CUDA
CUDA.allowscalar(true)
# size of array
dsize = 10_000_000
# index
idx_h = collect(1:Int64(dsize/2))
idx_d = CuArray(idx_h)
# datasets
dt1_h = ones(dsize)
dt2_d = CuArray(dt1_h)
dt3_d = CuArray(dt1_h)
dt4_d = CuArray(dt1_h);

```

## Functions:

1. cpu version

```julia
@views function t1!(idx_h, dt1_h)
    dt1_h[idx_h] .+= 1.0
    return nothing
end

```

2. gpu version: data and index are on GPU

```julia
@views function t2!(idx_d, dt2_d)
    dt2_d[idx_d] .+= 1.0
    return nothing
end

```

3. gpu version: index on cpu, data on GPU

```julia
@views function t3!(idx_h, dt3_d)
    dt3_d[idx_h] .+= 1.0
    return nothing
end

```

4. gpu version: kernel function

```julia
function kernel!(idx_d, dt4_d, sizeidx)
    ix = (blockIdx().x-1)*blockDim().x+threadIdx().x
    if ix≤sizeidx
        dt4_d[idx_d[ix]] += 1.0
    end
    return nothing
end
function t4!(idx_d, dt4_d)
    tds = 768
    bls = cld(length(idx_d), 768)
    CUDA.@sync begin
        @cuda threads=tds blocks=bls kernel!(idx_d, dt4_d, length(idx_d))
    end
    return nothing
end

```

## Results:

```julia-auto
@benchmark t1!($idx_h, $dt1_h)
@benchmark CUDA.@sync t2!($idx_d, $dt2_d)
@benchmark CUDA.@sync t3!($idx_h, $dt3_d)
@benchmark t4!($idx_d, $dt4_d)

```

| Item | Name | Time |
| --- | --- | --- |
| 1 | cpu version | 14.062 ms ± 2.516 ms |
| 2 | data and index are on GPU | 30.188 ms ± 10.241 ms |
| 3 | index on cpu, data on GPU | 18.524 ms ± 5.147 ms |
| 4 | kernel function | 163.619 μs ± 104.392 μs |

## Questions:

1. According to the results, putting the index on the gpu will be slower than `Item 3`. However, from my understanding, it seems that `Item2` should be a bit faster than `Item3` since they are both on the GPU.
2. I noticed that when I did the benchmark for `Item3`, GPU memory is almost fully occupied (24 GB). What leads to this phenomenon?
3. From the benchmark results it seems that kernel function is the most efficient way to execute, can we achieve similar performance with CuArray alone?

---

<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:** [November 23, 2022, 4:18pm UTC](https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600/2 "2022-11-23T16:18:10Z")

</div>

`view(dt4_d, idx_d) .+ 1` should essentially do the same as your custom kernel, by passing a SubArray to the broadcast kernel. So that involves much more complex functionality (SubArray and broadcast), but that shouldn’t warrant a 2 order of magnitude slowdown.

---

<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:** [November 23, 2022, 4:33pm UTC](https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600/3 "2022-11-23T16:33:25Z")

</div>

Oh wait, the problem is the memory copy that happens during `view` for the purpose of bounds checking. So if you do `@inbounds view(data, idx) .+= 1` that’s almost as fast as your custom kernel version.

---

<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:** [October 31, 2023, 10:29am UTC](https://discourse.julialang.org/t/correct-implementation-of-cuarrays-slicing-operations/90600/4 "2023-10-31T10:29:10Z")

</div>

This issue will be fixed by [Rework host indexing. by maleadt · Pull Request #499 · JuliaGPU/GPUArrays.jl · GitHub](https://github.com/JuliaGPU/GPUArrays.jl/pull/499). `@inbounds` will still be required to get to the level of performance of the custom kernel, but the bounds check is now significantly faster. See [view(data, idx) boundschecking is disproportionately expensive · Issue #1678 · JuliaGPU/CUDA.jl · GitHub](https://github.com/JuliaGPU/CUDA.jl/issues/1678#issuecomment-1786935526) for detailed timings.