# CUDA on slices

**URL:** <https://discourse.julialang.org/t/cuda-on-slices/86782>\
**Category:** General Usage\
**Tags:** cuda, broadcasting, cuarrays, tensors, mapslices\
**Created:** [September 5, 2022, 5:35am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782 "2022-09-05T05:35:28Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lincoln\_Hannah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lincoln_hannah/32/19198_2.png) [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)\
**Post date:** [September 5, 2022, 5:35am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/1 "2022-09-05T05:35:28Z")

</div>

Trying to apply a function to slices of a CuArray.  
Tried broadcasting, mapslices and various Tensor packages.  
Is there any way to do this?

```julia
using CUDA, Pipe, Interpolations, TensorCast

d = Uniform(1,30)
x = Float64.(1:30)
y = rand(d,30,500) |> eachcol
v = rand(d,1000)
gx, gy, gv = CuArray.([x,y,v])

f( y::Array{Float64}, v::Float64)::Float64 = linear_interpolation( x, y )(v)
gf(gy::CuArray{Float64}, gv::Float64)::Float64 = linear_interpolation( gx, gy )(gv)

@cast o[i,j] := f( y[i], v[j] ) #works
@cast go[i,j] := gf( gy[i], gv[j] ) #fails

```

---

<div class="post-metadata">

**Author:** ![Jarartur](https://avatars.discourse-cdn.com/v4/letter/j/b2d939/32.png) [@Jarartur](https://discourse.julialang.org/u/Jarartur)\
**Post date:** [September 5, 2022, 7:31am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/2 "2022-09-05T07:31:58Z")

</div>

> [@Lincoln\_Hannah](#):
>
> `@cast go[i,j] := gf( gy[i], gv[j] ) #fails`

What error does this give you? Is it scalar indexing? If so then it is a limitation of CUDA programming model as CUDA.jl does not currently automatically vectorize your functions which means that accessing a `CuArray[index]` will throw you an error if not explicitly allowed. You should be able to do `CuArray[1:100]` though.  
I am not familiar with TensorCast implementation and I don’t have a CUDA capable computer in this moment but you could look into that.  
If you want to do einsum on GPU I remember I used [Tullio](https://github.com/mcabbott/Tullio.jl) once and it worked.

---

<div class="post-metadata">

**Author:** ![Lincoln\_Hannah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lincoln_hannah/32/19198_2.png) [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)\
**Post date:** [September 5, 2022, 8:31am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/3 "2022-09-05T08:31:43Z")

</div>

I’ve paid for 4 more hours here. password Aa12345678  
[https://9jeqf.launch.juliahub.app/](https://9jeqf.launch.juliahub.app/)

```julia
using CUDA, Pipe, Interpolations, TensorCast, Distributions, Tullio

d = Uniform(1,30)
x = Float64.(1:30)
y = rand(d,30,500) |> eachcol |> collect
v = rand(d,1000)
gx, gy, gv = CuArray.([x,y,v])

f(y, v) = linear_interpolation( x, y )(v)
gf(gy, gv) = linear_interpolation( gx, gy )(gv)

@cast o[i,j] := f(y[i],v[j]) #works
@tullio o[i,j] := f(y[i],v[j]) #works

@cast go[i,j] := gf(gy[i],gv[j]) #error
@tullio go[i,j] := gf(gy[i],gv[j]) #error

```

@cast gives the error

```julia
ERROR: GPU broadcast resulted in non-concrete element type Any.
This probably means that the function you are broadcasting contains an error or type instability.

```

@tullio gives

```julia
Warning: Performing scalar indexing on task Task (runnable) @0x00007fc6df124740.
│ Invocation of getindex resulted in scalar indexing of a GPU array.
│ This is typically caused by calling an iterating implementation of a method.
│ Such implementations *do not* execute on the GPU, but very slowly on the CPU,
│ and therefore are only permitted from the REPL for prototyping purposes.
│ If you did intend to index this array, annotate the caller with @allowscalar.

```

---

<div class="post-metadata">

**Author:** ![Jarartur](https://avatars.discourse-cdn.com/v4/letter/j/b2d939/32.png) [@Jarartur](https://discourse.julialang.org/u/Jarartur)\
**Post date:** [September 5, 2022, 5:06pm UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/4 "2022-09-05T17:06:12Z")

</div>

For Tullio to work with GPU you need to import also:

```julia
using CUDA, CUDAKernels, KernelAbstractions

```

And from what I see it is some problem with your call to `linear_interpolation`.

```julia
julia> @tullio o[i,j] := linear_interpolation( gx, gy[i] )(gv[j])
Error: Reason: unsupported dynamic function invocation (call to linear_interpolation)

```

for more you can read [CUDA Troubleshooting](https://cuda.juliagpu.org/stable/development/troubleshooting/#Troubleshooting).  
If we dig deeper we can see it is some conversion function

```julia
Reason: unsupported dynamic function invocation (call to convert)

```

And if we inspect `gy` we find out that you didn’t really move y to CuArray

```julia
julia> typeof(gy[1])
SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}

```

_(see how it is `Matrix{Float64}` and not `CuArray{Float64}`, that’s because broadcasting works only one level deep.)_

This is the hardest part of GPU programming as GPUs don’t really like lists or Vectors of Vectors (CUDA.jl will throw:

```julia
ERROR: CuArray only supports element types that are stored inline

```

if you will try to do this.)

The “easiest” way to do this would be to write a custom kernel, or use einsum to account for two dimensions and store `gy` as a CuArray of size (30,500).

Sorry but here my knowledge ends unfortunately. Also I don’t know about your experience but keep in mind that Julia tries to compile a GPU kernel in contrast to i.e. Python frameworks like PyTorch or CuPy where the kernels are already compiled and you only use an api to access them, so in CUDA.jl you generally need to work on Arrays of numbers only.

---

<div class="post-metadata">

**Author:** ![Lincoln\_Hannah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lincoln_hannah/32/19198_2.png) [@Lincoln\_Hannah](https://discourse.julialang.org/u/Lincoln_Hannah)\
**Post date:** [September 6, 2022, 2:18am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/5 "2022-09-06T02:18:08Z")

</div>

Thanks Jaratur.  
I tried a few einsum packages. When they do work they tend to be slower than `mapslices` on a CPU.

I’ve reposted here [Mapslices very slow](https://discourse.julialang.org/t/mapslices-very-slow/86826)

I’d rather not write a custom Kernel. See if there’s any alternatives.

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [September 6, 2022, 2:45am UTC](https://discourse.julialang.org/t/cuda-on-slices/86782/6 "2022-09-06T02:45:33Z")

</div>

Take a look at this (open) [issue](https://github.com/JuliaMath/Interpolations.jl/issues/357#issuecomment-1154771303) and a related merged [PR](https://github.com/JuliaGPU/CUDA.jl/pull/460#issuecomment-701441464) for GPU interpolation. It’s immature, but better than trying to roll your own.