GPU#

Le calcul sur GPU permet de calculer plus rapidement certaines opérations mathématiques. Il est particulièrement bien adapté pour les opérations simples entre tableaux de grandes dimensions. Il consiste à déporter les calculs sur la carte graphique puis récupérer les résultats. Ces transferts de données ont un coût qu’il faudra prendre en compte lors de l’utilisation de ces co-processeurs.

On importe quelques packages pour nos exemples

using BenchmarkTools
using Random
using Test
using LinearAlgebra
using ForwardDiff
using ProgressMeter
using Plots

Le package principal qui permet d’utiliser est CUDA.jl, il permet d’utiliser les cartes graphiques de la marque NVIDIA.

using CUDA

CUDA.versioninfo()

CUDA toolchain: 
- runtime 12.9, artifact installation
- driver 575.57.8 for 12.9
- compiler 12.9

CUDA libraries: 
- cuBLAS: 12.9.1
- cuSPARSE: 12.5.10
- cuSOLVER: 11.7.5
- cuFFT: 11.4.1
- cuRAND: 10.3.10
- CUPTI: 2025.2.1 (API 12.9.1)
- NVML: 12.0.0+575.57.8

Julia packages: 
- CUDACore: 6.0.0
- GPUArrays: 11.5.3
- GPUCompiler: 1.9.1
- KernelAbstractions: 0.9.41
- CUDA_Driver_jll: 13.2.1+0
- CUDA_Compiler_jll: 0.4.3+0
- CUDA_Runtime_jll: 0.21.0+1

Toolchain:
- Julia: 1.12.6
- LLVM: 18.1.7

1 device:
  0: Tesla V100S-PCIE-32GB (sm_70, 31.724 GiB / 32.000 GiB available)

Lorsque vous installez CUDA.jl, une version du compilateur de NVIDIA sera également téléchargé sur votre poste. Il est possible d’utiliser une installation existante, c’est expliqué dans la documentation.

Cette première fonction permet de vérifier que le package est correctement installé et que vous disposez de la carte graphique adéquate.

CUDA.functional()

true

Vous pouvez également lister les matériels à votre disposition.

for device in CUDA.devices()
    @show capability(device)
    @show name(device)
end

capability(device) = v"7.0.0"
name(device) = "Tesla V100S-PCIE-32GB"

Pour tester votre installation vous pouvez également regarder le package GPUInspector.jl

Création de tableaux pour le GPU#

Allocations sur la carte graphique#

a = CuArray{Float32,2}(undef, 2, 2)

2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:
 0.0  0.0
 0.0  0.0

similar(a)

2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:
 0.0  0.0
 0.0  0.0

a = CuArray([1,2,3])

3-element CuArray{Int64, 1, CUDACore.DeviceMemory}:
 1
 2
 3

Transfert vers le CPU#

b est alloué sur le CPU, et un transfert de données est effectué. Deux possibilités sont offertes, transformer votre tableau en type Array ou utiliser la fonction collect.

b = Array(a)

3-element Vector{Int64}:
 1
 2
 3

collect(a)

3-element Vector{Int64}:
 1
 2
 3

Compatibilité avec les tableaux Julia#

CUDA.ones(2)

2-element CuArray{Float32, 1, CUDACore.DeviceMemory}:
 1.0
 1.0

a = CUDA.zeros(Float32, 2)

2-element CuArray{Float32, 1, CUDACore.DeviceMemory}:
 0.0
 0.0

a isa AbstractArray

true

CUDA.fill(42, (3,4))

3×4 CuArray{Int64, 2, CUDACore.DeviceMemory}:
 42  42  42  42
 42  42  42  42
 42  42  42  42

Tirages aléatoires#

CUDA.rand(2, 2)

2×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:
 0.936359  0.00342143
 0.274716  0.112303

CUDA.randn(2, 1)

2×1 CuArray{Float32, 2, CUDACore.DeviceMemory}:
  0.9938167
 -0.6869481

x = CUDA.CuArray(0:0.01:1.0)
nt = length(x)
y = 0.2 .+ 0.5 .* x + 0.1 .* CUDA.randn(nt);
scatter( Array(x), Array(y))
plot!( x -> 0.2 + 0.5x)
xlims!(0,1)
ylims!(0,1)

X = hcat(CUDA.ones(nt), x);

β = X'X \ X'y

2-element CuArray{Float64, 1, CUDACore.DeviceMemory}:
 0.1942306079021719
 0.5163624904090824

a = CuArray([1 2 3])

1×3 CuArray{Int64, 2, CUDACore.DeviceMemory}:
 1  2  3

view(a, 2:3)

2-element CuArray{Int64, 1, CUDACore.DeviceMemory}:
 2
 3

a = CuArray{Float64}([1 2 3])
b = CuArray{Float64}([4 5 6])

map(a) do x
    x + 1
end

1×3 CuArray{Float64, 2, CUDACore.DeviceMemory}:
 2.0  3.0  4.0

reduce(+, a)

6.0

accumulate(+, b; dims=2)

1×3 CuArray{Float64, 2, CUDACore.DeviceMemory}:
 4.0  9.0  15.0

findfirst(isequal(2), a)

CartesianIndex(1, 2)

a = CuArray([1 2 3])
b = CuArray([4 5 6])

map(a) do x
    x + 1
end

a .+ 2b

reduce(+, a)

accumulate(+, b; dims=2)

findfirst(isequal(2), a)

CartesianIndex(1, 2)

CUBLAS#

Les opérations entre tableaux alloués sur le GPU sont effectués sur le co-processeur.

a * b'

1×1 CuArray{Int64, 2, CUDACore.DeviceMemory}:
 32

CUSOLVER#

Certaines fonctions d’algèbre linéaire de LAPACK sont disponibles:

A =  cu([1.80   2.88   2.05  -0.89;
  5.25  -2.95  -0.95  -3.80;
  1.58  -2.69  -2.90  -1.04;
 -1.11  -0.66  -0.59   0.80])

B = cu([9.52  18.47;
 24.35   2.25;
  0.77 -13.28;
 -6.22  -6.21   ])

4×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:
  9.52   18.47
 24.35    2.25
  0.77  -13.28
 -6.22   -6.21

L, ipiv = CUDA.CUSOLVER.getrf!(A)
CUDA.CUSOLVER.getrs!('N', L, ipiv, B)

4×2 CuArray{Float32, 2, CUDACore.DeviceMemory}:
  1.0  3.0
 -1.0  2.0
  3.0  4.0
 -5.0  0.999996

CUFFT#

Pour faire des FFTs sur le GPU, il est nécessaire d’utiliser les plans pour allouer l’espace nécessaire sur la carte graphique.

fft = CUFFT.plan_fft(a) 
fft * a

1×3 CuArray{ComplexF64, 2, CUDACore.DeviceMemory}:
 6.0+0.0im  -1.5+0.866025im  -1.5-0.866025im

ifft = CUFFT.plan_ifft(a)
real(ifft * (fft * a))

1×3 CuArray{Float64, 2, CUDACore.DeviceMemory}:
 1.0  2.0  3.0

CUSPARSE#

Les formats destinés aux matrices creuses sont aussi disponibles.

A = SymTridiagonal(
                [ 4.0,  10.0,  29.0,  25.0 ,  5.0], [-2.0,  -6.0,  15.0 ,  8.0])

5×5 SymTridiagonal{Float64, Vector{Float64}}:
  4.0  -2.0    ⋅     ⋅    ⋅ 
 -2.0  10.0  -6.0    ⋅    ⋅ 
   ⋅   -6.0  29.0  15.0   ⋅ 
   ⋅     ⋅   15.0  25.0  8.0
   ⋅     ⋅     ⋅    8.0  5.0

using SparseArrays

A = CUDA.CUSPARSE.CuSparseMatrixCSR(sparse(A))

5×5 cuSPARSE.CuSparseMatrixCSR{Float64, Int32} with 13 stored entries:
  4.0  -2.0    ⋅     ⋅    ⋅ 
 -2.0  10.0  -6.0    ⋅    ⋅ 
   ⋅   -6.0  29.0  15.0   ⋅ 
   ⋅     ⋅   15.0  25.0  8.0
   ⋅     ⋅     ⋅    8.0  5.0

Les étapes pour porter votre code sur GPU#

  1. Développez votre application sur votre CPU avec les tabeaux de type Array

  2. Portez votre application sur GPU en utilisant le type CuArray

  3. Desactivez “scalar indexing” pour trouver les opérations incompatibles.

  4. Ecrivez vos propres kernels CUDA pour remplacer ces opérations.

Noyaux CUDA#

L’écriture directe de noyaux CUDA est possible, cependant:

  • les allocations sont interdites,

  • pas d’entrées-sorties donc pas d’affichage,

  • si votre code n’est pas typé correctement, le code compilé sera peu performant.

Programmer vos noyaux de manière incrémentale, en les gardant le plus simple possible et en vérifiant soigneusement que le résultat escompté est correct.

a = CUDA.zeros(1024)

function kernel(a)
    i = threadIdx().x
    a[i] += 1
    return
end

@cuda threads=length(a) kernel(a)

CUDACore.HostKernel for kernel(CuDeviceVector{Float32, 1})

a = CUDA.rand(Int, 1000)

1000-element CuArray{Int64, 1, CUDACore.DeviceMemory}:
  8024292473858245638
  8366351865986462909
  -478590473994982443
  6455011415532812414
  1971351246391272035
  5127364939284440565
  2630115178803638390
 -5644401552727578332
 -5140528744280786978
 -8499091271048549664
 -7694186532867547595
 -5255459634966745612
 -8333941916029631128
                    ⋮
  6957052522391132310
  7706424498257350667
  -511530605672646660
  1377069193052232014
   530995064234415733
  5972871715024742470
   873356819365331130
 -3572393798617098110
 -7632414263804863927
  5519996393447067089
   959190096788747155
 -2201510031648981104

norm(a)

1.715004771591657e20

@btime norm($a)

  83.401 μs (147 allocations: 4.17 KiB)

1.715004771591657e20

@btime norm($(Array(a)))

  1.704 μs (0 allocations: 0 bytes)

1.715004771591657e20

La fonction norm est bien plus rapide exécutée sur le CPU

CUDA.allowscalar(false)

a = CuArray(1:9_999_999);

@time a .+ reverse(a);

  0.003002 seconds (491 allocations: 12.594 KiB)

Pour effectuer cette dernière instruction, vous avez besoin de programmer deux noyaux. La macro @time n’est pas adéquate pour évaluer la performance car on a affaire à une opération de type “lazy”. C’est à dire que l’expression est programmée sur le GPU mais pas exécutée. Elle le sera lorsque vous transférerez le résultat vers le CPU. Vous pouvez utiliser @sync ou @time du package CUDA.

@time CUDA.@sync a .+ reverse(a);

  0.000578 seconds (223 allocations: 4.766 KiB)

CUDA.@time a .+ reverse(a);

  0.000847 seconds (196 CPU allocations: 4.344 KiB) (2 GPU allocations: 152.588 MiB, 44.69% memmgmt time)

@btime CUDA.@sync $a .+ reverse($a);

  424.485 μs (80 allocations: 2.39 KiB)

@btime CUDA.@sync $(Array(a)) .+ reverse($(Array(a)));

  92.583 ms (8 allocations: 152.59 MiB)

Aide au développement#

Vous avez quelques macros disponibles pour vous aidez à implémenter vos noyaux:

kernel() = (@cuprintln("foo"); return)

kernel (generic function with 2 methods)

@cuda kernel()

CUDACore.HostKernel for kernel()

kernel() = (@cuprintln("bar"); return)

kernel (generic function with 2 methods)

@cuda kernel()

foo

CUDACore.HostKernel for kernel()

Attention, exécuter plusieurs noyaux CUDA à la suite prend du temps car il y a un temps de lattence plus important que sur CPU.

a = CuArray(1:9_999_999)
c = similar(a)
c .= a .+ reverse(a);

bar

function vadd_reverse(c, a, b)
    i = threadIdx().x
    if i <= length(c)
        @inbounds c[i] = a[i] + reverse(b)[i]
    end
    return
end

vadd_reverse (generic function with 1 method)

Essayons de remplacer la fonction reverse

function vadd_reverse(c, a, b)
    
    i = threadIdx().x
    if i <= length(c)
        @inbounds c[i] = a[i] + b[end - i + 1]
    end
    return
end

vadd_reverse (generic function with 1 method)

Cela ne fonctionne pas car on n’itère pas sur un tableau alloué sur un GPU de la même manière qu’un tableau alloué sur CPU.

try
    @cuda threads = length(a) vadd_reverse(c, a, a)
catch e
    println(e)
end

Number of threads in x-dimension exceeds device limit (9999999 > 1024).

Les fonctions blockIdx et threadIdx sont là pour vous aider:

attribute(device(), CUDA.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)

1024

function vadd_reverse(c, a, b)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(c)
        @inbounds c[i] = a[i] + b[end - i + 1]
    end
    return
end

vadd_reverse (generic function with 1 method)

Le noyau construit est plus rapide que la fonction reverse initialement utilisée:

@btime CUDA.@sync @cuda threads=1024 blocks=length($a)÷1024+1 vadd_reverse($c, $a, $a)

  249.664 μs (19 allocations: 1.03 KiB)

CUDA.HostKernel for vadd_reverse(CuDeviceVector{Int64, 1}, CuDeviceVector{Int64, 1}, CuDeviceVector{Int64, 1})

@btime CUDA.@sync $a .+ reverse($a);

  422.743 μs (74 allocations: 3.98 KiB)