# Power Consumption Comparison of GPU Linear Solvers for Cellular Potts Model Simulations
## 2 Mathematical Background
In this case, the expression for \(K_{e}\) is given by  
\[K_{e}=\int_{\Omega_{e}}B^{T}DB\,d\Omega_{e} \tag{8}\]  
Set \(\Omega_{e}\) refers to the specific domain associated with element \(e\), while \(B\) is identified as the strain-displacement matrix in the case of quadrilateral elements with four nodes. Finally, the material matrix \(D\) is related to the plane stress conditions [12]. By denoting by \(K\) the global stiffness matrix, which is assembled by putting together the local stiffness matrices \(K_{e}\), the global mechanical behavior of the system is realized by solving a large set of linear equations, represented as  
\[Ku=f. \tag{9}\]  
Here, \(u\) represents the vector of all displacements of all nodes, and the vector  
\[f=(f_{P})_{P\in\Omega} \tag{10}\]  
stands for the forces \(f_{P}\) acting on all pixels \(P\). Following the description of [PERSON] and [PERSON] in [13], for each \(P\) the cumulative force \(f_{P}\) is the resultant vector of all forces applied by its neighboring pixels:  
\[f_{P}=\sum_{Q\in\operatorname{\mathbb{N}b}(P)}f_{P,Q}, \tag{11}\]  
with \(f_{P,Q}\) denoting the force applied to \(P\) by pixel \(Q\). The stiffness matrix \(K\) has to be assembled only once but at each time step it is employed in solving the linear system (9), which is a large-scale problem. In [9], the Jacobi iterative method, through the gesvdj routine from the cuSOLVER library, was employed as the solver. A study of its impact on power consumption was investigated in [10]. Here, a comparison, in terms of efficiency and energy saving, to three other linear solvers is provided: the LU decomposition, implemented through the getrf routine from cuSOLVER; the Gauss-Seidel method, implemented using the AmxG library [14]; and a custom-made version of the conjugate gradient method, which is particularly relevant since we are dealing with a symmetric and positive-definite matrix \(K\). Here, a short description of the CG method is provided. Starting from an initial guess \(u_{0}\), the CG iteratively computes approximation \(u_{k}\)\((k=1,2,\dots)\) by optimizing along mutually conjugate directions. Using the notation in (9), the process minimizes the value of the quadratic function  
\[F(u)=\frac{1}{2}u^{T}Ku-u^{T}f \tag{12}\]  
over the next orthogonal search directions. At each iteration, the approximation \(u_{k}\) is computed in the following way. First, a corresponding initial residual \(r_{0}=f-Ku_{0}\) to the initial guess \(u_{0}\) is computed. Then, after setting \(d_{0}=r_{0}\), for \(k=0,1,\dots\) (until convergence), the following steps are iteratively executed:  
(i) Compute:  
\[\alpha_{k}=\frac{r_{k}^{T}r_{k}}{d_{k}^{T}kd_{k}}\ ; \tag{13}\]  
(ii) Update approximation: \(u_{k+1}=u_{k}+\alpha_{k}d_{k}\);  
(iii) Update residual \(r_{k+1}=r_{k}-\alpha_{k}Kd_{k}\), and the direction  
\[H_{v}=\lambda\sum_{k}\left(\frac{V_{k}}{\bar{V}}-1\right)^{2},\quad V_{k}= \operatorname{vol}(C_{k}).\]0  
The above discussion is summarized in the following Algorithm 1, which is the base of the parallel implementations proposed in the next section.