Why is the ETA of simulate() on my CPU faster than my GPU?

[nathanjbeaumont]

Hi everyone,
I am simulating a pretty hefty Pulseq sequence with 1000 TRs and a spiral trajectory (1192276 time points). I made a phantom with 232050 spins and with my GPU (RTX 2070, 8GB RAM), with Nblocks=1,000 and precision="f32" I get an ETA of 3+ hours. However if I set "gpu" = false and Nthreads=12 my CPU (AMD Ryzen 5 3600 6-core processor) has an ETA of 15 minutes.

My intuition tells me the GPU sim should work faster than my CPU. Is there a setting I'm missing?

Thanks in advance

[cncastillo]

Which version of Koma are you using? It is true that the GPU will not always be faster than the CPU, specially if you have a powerful CPU and not so powerful GPU. In the newest versions we have a pretty optimized CPU method (>= 0.9, I think). The GPU performance benefits are seen when the number of spins increases, or by increasing the simulation block size. We may need a better method to control the later, now we specify the number of blocks, while a max block size might make more sense (number of time points per kernel call). You should definitely be in the regime where the GPUs shine, just by looking at the number of spins.

CPU: Bloch with CPU
Array-based: BlochSimple
Kernel-based: Bloch with GPU

Now, the time difference you are observing is kind of crazy (CPU 15 min vs GPU 3h). I have never seen anything like that, and I have simulated pretty long sequences. Can you give us a little more information about the simulation? or something I can test (like 1/10 of the sequence, as it should scale linearly with the sequence duration)?

Off-topic: really happy people are really pushing the simulations :) 232050 spins, and 1192276 points. Very cool! 😄

[nathanjbeaumont]

@cncastillo are you able to reproduce what I'm seeing on your machine?

[cncastillo]

Mmm not really. In Julia 1.10 with KomaMRI 0.9.1, I get:

CPU Intel(R) Xeon(R) Gold 6142 CPU @ 2.60GHz, 64 threads: (sim_params: "gpu" = false, "Nthreads" = Threads.nthreads())
- Results: 302.094539 seconds
GPU NVIDIA A100 80GB PCIe: (sim_params: "gpu" = true, "Nblocks" = 1_000)
- Results: 124.722699 seconds

I will try with 1.11 now.

[cncastillo]

Same results with Julia 1.11.6:

Julia 1.11.6
CPU Intel(R) Xeon(R) Gold 6142 CPU @ 2.60GHz, 64 threads: 
- Results: 301.086687 seconds
GPU NVIDIA A100 80GB PCIe:
- Results: 140.709110 seconds

I had to change the code to make it run. It was using a Python-like syntax for getting the size of an array x.array[i] instead of size(x, i).

using CUDA, FileIO, KomaMRI, MAT, Images
cd(@__DIR__)
t1_img = load("calimetrix_phantom_t1.png")
imgarr = channelview(t1_img)
imgarr = imgarr[1,:,:]
t1_img_float = Float64.(imgarr)
t1_img_float = repeat(t1_img_float, outer = (1,1,7))
T1 = vec(t1_img_float) * 1200e-3  # seconds

t2_img = load("calimetrix_phantom_t2.png")

imgarr_t2 = channelview(t2_img)
imgarr_t2 = imgarr_t2[1,:,:]
t2_img_float = Float64.(imgarr_t2)
t2_img_float = repeat(t2_img_float, outer = (1,1,7))
T2 = vec(t2_img_float) * 100e-3  # seconds

res = 125e-3 / size(t1_img_float, 2)  # get resolution in m from diameter
xs = ( -(size(t1_img_float, 1) - 1)/2:(size(t1_img_float, 1) - 1)/2 ) * res
ys = ( -(size(t1_img_float, 2) - 1)/2:(size(t1_img_float, 2) - 1)/2 ) * res
zs = ( -(size(t1_img_float, 3) - 1)/2:(size(t1_img_float, 3) - 1)/2 ) * 1e-3

x_coords = float(vec([x for x=xs, y=ys, z=zs]))
y_coords = float(vec([y for x=xs, y=ys, z=zs]))
z_coords = float(vec([z for x=xs, y=ys, z=zs]))
ρ = ones(size(z_coords, 1))

obj = Phantom(name="spin", x=x_coords, y=y_coords, z=z_coords, ρ=ρ, T1=T1, T2=T2, T2s=T2)

seq_name = "mrf_spiralG_30cm_300SR_70G"
seq = read_seq("$(seq_name).seq")

##
sys = Scanner(B0=1.5, Smax=333, Gmax=80e-3)
sim_params = KomaMRICore.default_sim_params()
sim_params["Nthreads"] = Threads.nthreads()
sim_params["Nblocks"] = 1_000
sim_params["precision"] = "f32"
sim_params["gpu"] = true
sim_params["return_type"] = "mat"
raw_mat = simulate(obj, seq, sys; sim_params)
# CUDA.reclaim() # <-- This is not necessary 
filename = "koma_sim_$(seq_name)"
matwrite("$(filename).mat", Dict(filename => raw_mat))

[nathanjbeaumont]

OK on my system the GPU ETA is 3 hours at NBlocks = 1_000, and plateaus down to 1 hour at NBlocks = 8_000. I waited a few minutes each try to confirm the ETA was accurate.

[cncastillo]

Can you wait for one of them to finish? To see if it actually takes 3h, based on my experience it will take a lot less. I want to see if it is a problem of the simulation speed or the ETA in your system.

[cncastillo]

The problem might be that the ETA, for the first iterations it considers precompilation, and then it slowly gets corrected. But the simulation time is actually much less with the GPU.