# How to properly pass structs into GPU? (MWE included)

**URL:** <https://discourse.julialang.org/t/how-to-properly-pass-structs-into-gpu-mwe-included/93727>\
**Category:** GPU\
**Created:** [January 29, 2023, 9:53am UTC](https://discourse.julialang.org/t/how-to-properly-pass-structs-into-gpu-mwe-included/93727 "2023-01-29T09:53:45Z")\
**Posts on this page:** 1\
**Showing post:** 7

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [January 29, 2023, 8:23pm UTC](https://discourse.julialang.org/t/how-to-properly-pass-structs-into-gpu-mwe-included/93727/7 "2023-01-29T20:23:50Z")

</div>

Thank God, I was finally able to do it exactly how I wanted! And this answers my 5th question:

_5. Is the way shown here the intended and “state of the art” approach for using structs in GPU programming - if not, what is? I looked at the interpolate example in the docs, and I struggle with grasping it._

Short answer, no. The code below is in my opinion much better:

```julia
# Packages
using CUDA
using Parameters
using StaticArrays
using Adapt
using Setfield
using LinearAlgebra
using BenchmarkTools

# Constants
const T = Float32
const ST = 3
const H = T(0.5)
const αD = T((7/(4*π*H^2)))
const ϵ = T(1e-6)
const NP = 10^4
const NL = 10^7

# Main Struct for Gradient Values
@with_kw struct ∇ᵢWᵢⱼStruct{T,ST}
    NL::Base.RefValue{Int}

    # Input for calculation
    xᵢⱼ ::AbstractVector{SVector{ST,T}} = zeros(SVector{ST,T},NL[])
    xᵢⱼ² ::AbstractVector{T} = similar(xᵢⱼ,T,NL[])
    dᵢⱼ ::AbstractVector{T} = similar(xᵢⱼ,T,NL[])
    qᵢⱼ ::AbstractVector{T} = similar(xᵢⱼ,T,NL[])
    ∇ᵢWᵢⱼ ::AbstractVector{SVector{ST,T}} = similar(xᵢⱼ,NL[])
end

# Helper struct to convert onto GPU
struct ∇ᵢWᵢⱼDeviceStruct{ST,T}
    xᵢⱼ ::CuDeviceVector{SVector{ST,T}, 1}
    xᵢⱼ² ::CuDeviceVector{T, 1}             
    dᵢⱼ ::CuDeviceVector{T, 1}             
    qᵢⱼ ::CuDeviceVector{T, 1}             
    ∇ᵢWᵢⱼ ::CuDeviceVector{SVector{ST,T}, 1} 
end

# Adapt Structure for Conversion to GPU
Adapt.adapt_structure(to, s::∇ᵢWᵢⱼStruct) = ∇ᵢWᵢⱼDeviceStruct(
                                                        adapt(to, s.xᵢⱼ ),
                                                        adapt(to, s.xᵢⱼ² ),
                                                        adapt(to, s.dᵢⱼ ),
                                                        adapt(to, s.qᵢⱼ ),
                                                        adapt(to, s.∇ᵢWᵢⱼ),
                                                        )

# Simple Grad Version
function GradVersion0(V::∇ᵢWᵢⱼStruct)
    V.xᵢⱼ² .= dot.(V.xᵢⱼ,V.xᵢⱼ)
    V.dᵢⱼ .= sqrt.(V.xᵢⱼ²)
    V.qᵢⱼ .= V.dᵢⱼ ./ H

    @. V.∇ᵢWᵢⱼ = αD *5 * (V.qᵢⱼ- 2)^3 * V.qᵢⱼ / (8H * (V.qᵢⱼ * H + ϵ)) * V.xᵢⱼ

    return nothing
end

### --------------------------- Writing Our Own Custom Kernel ---------------------------

# Actual Calculation of Gradient in Kernel
function GPU_GRAD(V)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x
    for i = index:stride:length(V.xᵢⱼ)
        xdot = dot(V.xᵢⱼ[i],V.xᵢⱼ[i])
        @inbounds V.xᵢⱼ²[i] += xdot
        dist = sqrt(xdot)
        @inbounds V.dᵢⱼ[i] += dist
        q = dist/H
        @inbounds V.qᵢⱼ[i] += q
        wg = αD * 5 * (q - 2)^3 * q / (8H * (q * H + ϵ)) * V.xᵢⱼ[i]
        @inbounds V.∇ᵢWᵢⱼ[i] += wg
    end
    return nothing
end

# Actual Function to Call for Kernel Gradient Calc
function GradVersion0_Kernel(V)
    kernel = @cuda launch=false GPU_GRAD(V)
    config = launch_configuration(kernel.fun)
    threads = min(length(V.xᵢⱼ), config.threads)
    blocks = cld(length(V.xᵢⱼ), threads)

    CUDA.@sync begin
        kernel(V; threads, blocks)
    end
end

### --------------------------- Benchmarking ---------------------------
IniX = rand(SVector{ST,T},NL)
V_CPU = ∇ᵢWᵢⱼStruct{T,ST}(NL=Ref(NL),xᵢⱼ=IniX)
V_GPU_Simple = ∇ᵢWᵢⱼStruct{T,ST}(NL=Ref(NL),xᵢⱼ=CuArray(IniX))
V_GPU_Kernel = ∇ᵢWᵢⱼStruct{T,ST}(NL=Ref(NL),xᵢⱼ=CuArray(IniX))

@CUDA.time GradVersion0(V_CPU)
@CUDA.time GradVersion0(V_GPU_Simple)
@CUDA.time GradVersion0_Kernel(V_GPU_Kernel)

if isapprox(V_CPU.∇ᵢWᵢⱼ,Array(V_GPU_Simple.∇ᵢWᵢⱼ), rtol=1e-6) & isapprox(V_CPU.∇ᵢWᵢⱼ,Array(V_GPU_Kernel.∇ᵢWᵢⱼ), rtol=1e-6)
    println("All three approaches produce the same result with rtol=1e-6.")
end

V_CPU_Bench = @benchmark GradVersion0($V_CPU)
# I could not get sensible timings other than including @CUDA.sync here too, and not only in GradVersion0_Kernel
V_GPU_Simple_Bench = @benchmark @CUDA.sync GradVersion0($V_GPU_Simple)
V_GPU_Kernel_Bench = @benchmark @CUDA.sync GradVersion0_Kernel($V_GPU_Kernel)

println("V_CPU_Bench")
display(V_CPU_Bench)
println("V_GPU_Simple_Bench")
display(V_GPU_Simple_Bench)
println("V_GPU_Kernel_Bench")
display(V_GPU_Kernel_Bench)

```

The code now consists of one immutable struct `∇ᵢWᵢⱼStruct` and its helper struct to turn it to the GPU using Adapt, `∇ᵢWᵢⱼDeviceStruct`. This is done through the “adapt\_structure” functionality highlighted in the link I referenced in my first post.

The really important part of `∇ᵢWᵢⱼStruct` is that it is now basedd on `AbstractArray` instead of being limited to CPU/GPU specifically - this allows for allocating the structure on both CPU and GPU and run the exact same code on it, so that is really great! - Also what the benchmark example does, for 10\*10^6 elements we see:

```julia
  0.104691 seconds (26 CPU allocations: 1.172 KiB)
  0.015593 seconds (225 CPU allocations: 17.391 KiB)
  0.013344 seconds (80 CPU allocations: 4.875 KiB)
All three approaches produce the same result with rtol=1e-6.
V_CPU_Bench
BenchmarkTools.Trial: 69 samples with 1 evaluation.
 Range (min … max): 69.430 ms … 105.620 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 72.026 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 73.314 ms ± 4.907 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

     ▂█▁ ▁ ▂  
  ▅▁██████▁▆█▅█▃▃▃▅▅▆▁▃▁▅▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▁
  69.4 ms Histogram: frequency by time 89.4 ms <

 Memory estimate: 1.17 KiB, allocs estimate: 26.
V_GPU_Simple_Bench
BenchmarkTools.Trial: 1863 samples with 1 evaluation.
 Range (min … max): 1.934 ms … 56.218 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 2.164 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 2.668 ms ± 2.434 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

  ██▅▄▄▄▂
  █████████▇▆▅▅▅▆▄▅▅▆▅▅▄▄▄▄▄▅▅▁▁▅▄▄▄▄▄▅▁▁▁▄▄▄▄▄▁▁▁▁▄▄▁▁▁▁▁▁▄ █
  1.93 ms Histogram: log(frequency) by time 12.5 ms <

 Memory estimate: 17.27 KiB, allocs estimate: 222.
V_GPU_Kernel_Bench
BenchmarkTools.Trial: 2343 samples with 1 evaluation.
 Range (min … max): 1.732 ms … 83.638 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 1.969 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 2.123 ms ± 1.767 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

         ▅██▄     
  ▂▃▄▄▆▇█████▇▄▄▃▃▃▃▃▃▃▃▃▃▃▂▂▃▃▃▃▃▃▃▃▃▂▂▃▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂ ▃
  1.73 ms Histogram: frequency by time 3.21 ms <

 Memory estimate: 4.84 KiB, allocs estimate: 79.

```

So this is extremely cool, being able to do both CPU/GPU with no hassle!

On top of this now the GPU functions accept the full struct with no problem due to:

```julia
isbits(cudaconvert(V_GPU_Kernel)) # true

```

So that is really awesome, so happy.

Thank you to @mcabbott for helping me understanding how to do this with an immutable struct.

Thank you to @lodygaw for the initial question and solution for how to do this, in a way I could understand.

So now I can get great performance, easy code to structure and organise - Julia is so cool right now 🙂

---

_[View the full topic](https://discourse.julialang.org/t/how-to-properly-pass-structs-into-gpu-mwe-included/93727)._