Add GPU reconstruction example to documentation

[nHackel]

[alexjaffray]

For anyone looking for how to perform the recon on the GPU.

To use the GPU functionality, you must add and use a GPU array library package of your choice (for example CUDA, ROCm, Metal etc.) in your environment i.e

Pkg.add(CUDA)

and recompile. If there are no unresolvable dependencies, then the MRIRecoGPUArraysExt module will be added for use. You don't have to explicitly add a use statement for it.

Then, in the parameter dictionary (for the current example let's call it reco_params) passed to the reconstruction call, add the parameter

reco_params[:arrayType] = CuArray

to use GPU and

reco_params[:arrayType] = Array

to use CPU (this is the default if the parameter is not given explicitly).

Do not try to specify the type beyond above, for example setting it as CuArray{ComplexF32,...} where that's the explicit type of your array does not work.

It's quite straightforward, hopefully this gives people a start before the docs are ready. I understand the interfaces and backends are currently undergoing changes (huge thanks to @nHackel and the team for putting effort into this) so it doesn't make sense to change the examples until things are finalized.

I can confirm, depending on matrix size and hardware, that GPU computation can be significantly faster than CPU.

For a test 96 x 96 x 72 3D reconstruction using 30 iterations and CGNR with L2 Regularization, with 16 channels, computation on CPU (32 threads @ 2.2 GHz) takes 3.5 seconds and computation on GPU (RTX A4500) takes 0.55 seconds.

[nHackel]

@alexjaffray thank you for the writeup, that's exactly the correct way to utilize GPU acceleration with MRIReco!

One further parameter which might be interesting to play around with is the scheduler argument which was added very recently. This allows you to change if seperate inverse problems of a reconstruction are solved in parallel or sequentially. The default for GPU based reconstructions is sequential, but depending on the size of the problem further improvements might be achievable with "further parallelism.

we use schedulers from OhMyThreads.jl, to change the scheduler we can either provide a symbol or a scheduler:

params[:scheduler] = DynamicScheduler() # or :dynamic for parallel execution
params[:scheduler] = SerialScheduelr() # or :serial for serial execution

I haven't played around much with different reconstruction sizes, so I don't have a good feeling when its helpful to change to dynamic problems. I'd image it would help if the problems are large enough to benefit from GPU acceleration and small enough for the GPU to handle multiple "instances" of the problem without fighting over resources. This sweet spot then of course depends on the GPU in question.

If a system has multiple GPUs to chose from, you can do that with CUDA.device!(...), there currently is no MRIReco specific thing to specify this for different reconstructions or distribute the inverse problems across GPUs (yet).

The interface and backend are actually fairly stable at the moment and most GPU related changes I'm thinking about are hopefully transparent to the user. So we could update the docs, I just didnt find the time yet 😅

One larger change which would require a docs rewrite is using AbstractImageReconstruction.jl as the interface for our reconstruction algorithms. This would allow us to make the algorithms more "modular" and extendable, but at the moment we aren't limited yet by the current interface like we were in MPIReco.jl

[alexjaffray]

All good! Thanks for the detail about the scheduler @nHackel ! I'll play around some more with the scheduler and see how it changes the reconstruction times when I have different matrix sizes and different dimensions.