# Power Consumption Comparison of GPU Linear Solvers for Cellular Potts Model Simulations
## 3 Parallel Approach  
In this section, the development of a parallel variant of Algorithm 1 is presented. The initial algorithm was implemented as sequential code in C. To address the complexity of the problem more effectively, the main components of this original CPU-based implementation were parallelized. While the sequential implementation delivered promising results, the related computational complexity and problem dimensions required an efficient solution, thus leading to the design of a parallel algorithm that mimics the fundamental operations of the sequential approach. Our implementation engages GPU technology to comprehensively demand its computational power across all components of the algorithm. This includes the initialization of cells, the construction of the lattice grid, and the implementation of the Metropolis algorithm. From the model description in Section 2 and Algorithm 1, it is clear that the most computationally intensive process, due to its frequent repetition, is to find the solution of the linear system (9). To address this, three different solvers are proposed, LU decomposition, Gauss-Seidel, and conjugate gradient (CG), in addition to the Jacobi method described in [9]. Given the properties of the matrix \(K\), CG is the most suitable. Hence, the focus is on the parallel implementation of the CG method, which is not included in the cuSOLVER library routines [15]. For precise control over how the algorithm is executed, all CUDA kernels are manually developed instead of using third party libraries. This choice enables the adaptation of the kernels exactly to specifications by optimizing performance and allowing for custom solutions. Following the structure outlined in Algorithm 1, a parallelization of the entire process using essential GPU-native API routines arising from CUDA [16] was implemented. This enables the user to optimize each stage for high performance. In order to discuss the parallel implementation, the main operations are split into three major steps, each designed to work as follows:  
**STEP 0** : Starting environment at time \(\tau=0\). In this step, the initialization of the lattice grid structure and cells using a thread-based parallelism strategy is addressed. This approach deploys pools of threads to manage the initialization of different cells segments. More precisely, a specific segment of a cell is assigned to one thread by maximizing the effective processing power of the GPU. This ensures efficient usage of GPU capabilities for priming the simulation environment. Adopting this approach, data transfers between the CPU and GPU are minimized, thereby reducing overhead typically associated with such transfers. Each thread is tasked with initializing a cell segment and assigning a unique identifier that defines its starting state. Additionally, these threads play a key role in building the lattice grid, linking each grid point to its respective cell. In this way, the initialization times are notably shortened. **STEP 1: FEM priming. ** The creation of matrix \(K\) through the FEM discretization process for strain is a one-time operation. To accomplish this efficiently, a dedicated CUDA kernel for constructing the global stiffness matrix was developed. This involves incorporating contributions from every mesh element. This process proceeds concurrently, with individual threads assigned to handle distinct segments of elements. This specific parallelization method guarantees extensive parallelism. **STEP 2: CPM iteration. ** This step collects specific custom CUDA kernels for each phase of the Metropolis algorithm: to update the force vectors and to assemble them, as in (10); to solve the linear system (9); and to calculate the change in energy, as in (7). In more detail:  
* Updating the force vector \(f\) is essential in simulations, as it reflects the active forces affecting each pixel and its neighboring pixels.