## 3 Parallel Programming and Implementations of GRICS
### Parallel Implementations of GRICS  
In this work, three different implementations were considered.  
* _OpenMP implementation_ In the pure OpenMP version of GRICS, the \(N_{c}\) coils are distributed over OpenMP threads by adding pragmas before the loop over coils, line 6 in Algorithm 2. The code is linked to threaded and serial FFTW-3.3.3 package [24], and to the Blitz-0.9 library [30]. There is no nested OpenMP parallelism in GRICS so far. The version of OpenMP is 4.5.2.
* _Hybrid MPI OpenMP implementation_ In the hybrid MPI and OpenMP version of the GRICS, the \(N_{c}\) coils are distributed over OpenMP threads, and the \(N_{t}\) motion states are distributed using the MPI, line 3 in Algorithm 2. Each MPI process represents an independent computational entity, operating on its allocated subset of motion states \(N_{t}\). The MPI guarantees inter-process communication and synchronizes computations among distributed processes or ranks. This distributed approach enables efficient workload allocation across the compute node. The OpenMP distribution is created explicitly by adding specific pragmas before the loop over coils. The code is linked to threaded and serial FFTW-3.3.3 package [24], and to the Blitz-0.9 library [30]. The version of the MPI is 4.0.3.
* _Preliminary MATLAB-GPU implementation of the reconstruction kernel_ In order to draw a comparison with CPU performance, we suggested a preliminary MATLAB-GPU implementation of the reconstruction kernel--the application of the operator \(E^{H}E\), step 6 of Algorithm 1, with iterative calls to Algorithm 2. For this purpose, the MATLAB _gpuArray_ object is used; it enables the storage of arrays on GPU memory, which allows the programmer to run a kernel on GPU with minimal software modifications. In particular, operators \(C,S,F\) of Algorithm 2 were declared as _gpuArray_. The spatial transformation operator \(T_{t}\) was the only operator that needed to be transferred from CPU to GPU memory at the beginning of each outer loop since their storage for all motion states \(t=1,2,\ldots,N_{t}\) on the GPU device memory resulted in an out-of-memory error on Quadro RTX 5000, which was not the case on A100). In MATLAB, the _gpuArray_ object supports CUDA, allowing MATLAB to efficiently harness the parallel processing capabilities of GPUs. By using CUDA, MATLAB can offload computationally intensive tasks to the GPU, resulting in speedups for certain types of calculations, particularly those involving large datasets. Additionally, within the framework of the MATLAB Parallel Computing Toolbox, the MAGMA library [31] is used to enhance the performance and accuracy of linear algebra functions designed for _gpuArray_ objects. Furthermore, MATLAB's execution model incorporates Just-In-Time (JIT) compilation, where code is compiled dynamically during runtime, optimizing performance by converting high-level MATLAB code into machine code. Although this approach is not as efficient as a full native GPU implementation, this JIT compilation process allows MATLAB to apply native GPU libraries to data stored in GPU memory, eliminating unnecessary data transfers between CPU memory and GPU memory. Additionally, we provided as Supplementary Materials a benchmarking of 3D FFT, which is at the core of all iterative MRI reconstruction programs. This provides an estimate of the upper bound of performance that could be achieved if the whole problem could fit in GPU memory.