# Power Consumption Comparison of GPU Linear Solvers for Cellular Potts Model Simulations
## 3 Parallel Approach
Initially, this process followed a sequential approach: beginning with the analysis of individual pixels, proceeding to identify their neighboring pixels, and concluding with the computation of specific pixel forces. While functional, this sequential method can be computationally intensive, especially in simulations involving large-scale systems. In contrast, the parallel implementation employs CUDA kernels to generate a matrix where each row corresponds to a unique pixel. This allows for simultaneous computation of forces across all pixels, eliminatingthe need for sequential pixel-by-pixel force calculations. By leveraging the parallel processing capabilities of GPU architectures, this approach maximizes efficiency through concurrent computations. Once forces for each pixel are asynchronously computed, they are then assigned to their respective slots in the force vector \(f\). This asynchronous process enhances performance by overlapping computational and memory operations, leveraging a critical advantage of the CUDA environment. * In order to solve the linear system presented in Equation (9), well-known solvers available in [14,15] (such as Jacobi, LU, and Gauss-Seidel) were considered. The choice of the conjugate gradient (CG) method was driven by the aim to investigate computational time and power consumption reductions. The CG method is inherently more efficient for the kind of symmetric and positive-definite matrices typically encountered in the proposed scheme, offering faster convergence rates compared to the Jacobi method. Moreover, an advantage of CG implementation is its fully GPU-centric design. Unlike the other methods, where employing routines from the cuSOLVER and AmgX libraries requires intermittent CPU-GPU data transfers, this approach fully operates on the GPU. This avoidance of data transfer not only notably reduces the overhead associated with these transfers but also speeds up the overall computational process. In the proposed parallel GPU implementation of the CG method, key operations such as vector addition (\(u_{k+1}=u_{k}+\alpha_{k}\cdot d_{k}\)) and subtraction (\(r_{k+1}=r_{k}-\alpha_{k}\cdot K\cdot d_{k}\)), as well as matrix-vector multiplication (\(K\cdot d_{k}\) and \(K\cdot u_{0}\)), are ad hoc designed for GPU execution. Basic algebraic operations are performed according to _dynamic parallelism_, equipped by most recent GPUs, which allows for nesting of CUDA kernels by enabling the invocation of other kernels from already executing threads. These threads handle the following operations: Vector computations are processed in a point-wise way, with each element processed independently by a dedicated thread by speeding up computation; the matrix-vector product is managed by distributing matrix rows across multiple threads, each one performing the dot product for its assigned row. In this way, the computations for the sparsity and structure of matrix \(K\) are optimized. In a similar way, the computation of the step size (\(\alpha_{k}\)) and direction vector (\(d_{k+1}\)) update also benefit from parallelization, using reduction patterns for dot products and point-wise operations for scalar-vector multiplications and additions. These CUDA kernel-based operations ensure coalesced memory access and use shared memory to minimize latency, modeling a fully based GPU approach. In addition, the CUDA threads' configuration (blocks\(\times\) threads) is defined in an _adaptive_ way according to the starting dimension of the \(\Omega\) domain. * The computation of energy changes is pivotal in the Metropolis algorithm. Three distinct CUDA kernels were developed to concurrently calculate the \(\Delta H_{ad}\), \(\Delta H_{v}\), and \(\Delta H_{s}\) components, employing an asynchronous parallel approach. Computing \(\Delta H_{ad}\) and \(\Delta H_{v}\) is straightforward and allows for rapid processing, thus these tasks are assigned to a single thread for optimal efficiency.