# Combine CartesianIndices for effective CUDA kernels

**URL:** <https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958>\
**Category:** GPU\
**Tags:** gpu, indexing\
**Created:** [May 11, 2021, 2:51pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958 "2021-05-11T14:51:25Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [May 11, 2021, 2:51pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958/1 "2021-05-11T14:51:25Z")

</div>

Suppose, I have a multidimensional array and I want to perform some operation along one of its axes. A sample generic code can look like this:

```julia
s = (1, 2, 3, 4, 5) # shape of the array
dim = 3 # dimension, along which I want to perform the operation

A = zeros(s)

spre = s[1:dim-1]
sdim = s[dim]
spost = s[dim+1:end]

Cpre = CartesianIndices(spre)
Cpost = CartesianIndices(spost)

for Ipost in Cpost
for Ipre in Cpre
    for i=1:sdim
        A[Ipre, i, Ipost] = i
    end
end
end

```

Next, I would like to launch these `for` loops on GPU using CUDA. To do it in a most efficient way I have to merge the two outer loops into one with an index through which I can iterate with strides. For that reason I need to merge the two Cartesian indices `Cpre` and `Cpost` into one set of indices. On CPU I came to the following solution:

```julia
Cprepost = CartesianIndices((spre..., spost...))

for I in Cprepost
    Ipre = Tuple(I)[1:dim-1]
    Ipost = Tuple(I)[dim:end]
    for i=1:sdim
        A[Ipre..., i, Ipost...] = i
    end
end

```

By analogy, I would expect that the kernel for GPU will look like

```julia
using CUDA

function kernel(A, dim, Cprepost, sdim)
    id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    for I in Cprepost[id:stride:end]
        Ipre = Tuple(I)[1:dim-1]
        Ipost = Tuple(I)[dim:end]
        for i=1:sdim
            A[Ipre..., i, Ipost...] = i
        end
    end

    return nothing
end

Agpu = CUDA.zeros(s)

@cuda threads=prod(spre)*prod(spost) kernel(Agpu, dim, Cprepost, sdim)

```

However, this code causes “unsupported call through a literal pointer” error.  
The only kernel which so far works for me is the following one:

```julia
function kernel(A, dim, Cprepost, sdim)
    id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    dim = 3

    for k=id:stride:length(Cprepost)
        Ipre = Tuple(Cprepost[k])[1:dim-1]
        Ipost = Tuple(Cprepost[k])[dim:end]
        for i=1:sdim
            A[Ipre..., i, Ipost...] = i
        end
    end
    return nothing
end

@cuda threads=prod(spre)*prod(spost) kernel(Agpu, dim, Cprepost, sdim)

```

However, here I have to explicitly define the `dim` variable within the kernel, because otherwise I obtain the same “unsupported call through a literal pointer” error.

Can you please help me to write the corresponding CUDA kernel.  
Thank you.

---

<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:** [May 11, 2021, 3:02pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958/2 "2021-05-11T15:02:34Z")

</div>

Have a look at CUDA.jl’s `mapreducedim!` kernel, it does things like that.

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [May 11, 2021, 3:02pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958/3 "2021-05-11T15:02:48Z")

</div>

> [@fedoroff](#):
>
> `Cprepost[id:stride:end]`

This creates a CPU array which can’t be dynamically allocated on the GPU. Try `@view Cprepost[id:stride:end] ` instead.

> [@fedoroff](#):
>
> However, here I have to explicitly define the `dim` variable within the kernel, because otherwise I obtain the same “unsupported call through a literal pointer” error.

You can use a [`Val` type](https://docs.julialang.org/en/v1/manual/types/#%22Value-types%22) to specialize on `dim`, which should make this type stable.

---

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [May 11, 2021, 3:36pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958/4 "2021-05-11T15:36:16Z")

</div>

It seems none of the approaches works.

P.S. I am reading CUDA sources now

---

<div class="post-metadata">

**Author:** ![fedoroff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fedoroff/32/53209_2.png) [@fedoroff](https://discourse.julialang.org/u/fedoroff)\
**Post date:** [May 11, 2021, 5:30pm UTC](https://discourse.julialang.org/t/combine-cartesianindices-for-effective-cuda-kernels/60958/5 "2021-05-11T17:30:46Z")

</div>

Well, it seems I found the solution:

```julia
using CUDA

s = (1, 2, 3, 4, 5)
dim = 3

spre = s[1:dim-1]
sdim = s[dim]
spost = s[dim+1:end]

Cpre = CartesianIndices(spre)
Cpost = CartesianIndices(spost)
C = CartesianIndices((length(Cpre), length(Cpost)))

function kernel(A, Cpre, Cpost, C, sdim)
    id = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for k=id:stride:length(C)
        i1 = C[k][1]
        i2 = C[k][2]
        Ipre = Cpre[i1]
        Ipost = Cpost[i2]
        for i=1:sdim
            A[Ipre, i, Ipost] = i
        end
    end
    return nothing
end

Agpu = CUDA.zeros(s)

@cuda threads=length(C) kernel(Agpu, Cpre, Cpost, C, sdim)

```

The trick was to create an additional set of Cartesian indices `C` which maps the linear index into the corresponding indices of `Cpre` and `Cpost`.