# Multi-GPU Acceleration for Finite Element Analysis in Structural Mechanics
## 3 Aggregation Algebraic Multigrid (AMG)
The resulting operators should satisfy \(B_{l}=\tilde{P}_{l}B_{l+1}\) and \(\tilde{P}_{l}^{T}\tilde{P}_{l}=I\), which implies that near null-space candidates lie in the rank of \(\tilde{P}_{l}\) and that columns of \(\tilde{P}_{l}\) are orthonormal [56]. We can use a smoother \(\tilde{S_{l}}\) to obtain a more robust prolongation operator \(P_{l}\) by smoothing the tentative prolongation matrix \(\tilde{P}_{l}\). The weighted Jacobi smoother is a common choice for this purpose. Finally, we perform the sparse Galerkin product as two matrix-matrix products. We use the AMG method as a preconditioner of a conjugate gradient iterative solver. Such an iterative solver is appropriate for linear systems of equations whose coefficient matrix is both positive-definite and symmetric as in the case of elasticity problems. Algorithm 2 shows the pseudo-code of the V-cycle of AMG for preconditioning the iterative solver. The algorithm aims to approximate the solution of (1) given the residual of the previous estimation of the iterative solver. The procedure consists of \(\mu_{1}\) smoothing operations to the approximate solution \(s_{l}\) at level \(l\) and the computation of the residual \(r_{l}\) for the relaxed approximate solution \(s_{l}\). This residual is then restricted to the coarse grid and solved if we reach the last \(L\) level. When we solve the approximate solution \(s_{l}\) in the coarsest grid, it is then prolonged to the finer grid by applying \(\mu_{2}\) smoothing operations to the approximate solution. Algorithm 3 details the conjugate gradient algorithm using the v-cycle of AMG as a preconditioner. This iterative solver only requires the maximum number of iterations \(\max_{\text{iter}}\), the tolerance \(<\mathit{tol}\), \(\mathrm{tol}_{\text{abs}}>\), the coefficient matrix \(\mathbf{K}\), the right-hand side \(\mathbf{f}\), and the initial guess for initializing the iterative procedure. The recursive solver provides an approximate solution \(\mathbf{u}\) of (1) with a residual after \(it\) iterations. ```
Input:\(\mathbf{n}_{1},\mathbf{s}_{1},1\) Output:\(\mathbf{s}_{1}\) // Pre-smooth/-relaxation: smooth\(\mu_{1}\) times on \(\mathbf{K}_{1}\)\(\mathbf{s}_{1}=\mathbf{r}_{1}\) \(\mathbf{s}_{1}\leftarrow\)\(\mathbf{Smooth}(\mathbf{r}_{1}\),\(\mathbf{K}_{1}\),\(\mathbf{s}_{1}\),\(\mu_{1})\) ; \(\mathbf{r}_{1}\leftarrow\mathbf{r}_{1}\) - \(\mathbf{K}_{1}\) - \(\mathbf{K}_{1}\) ; // Compute residual \(\mathbf{r}_{1+1}\gets R_{I}\) ; // Restrict residual to coarse grid // Coarsest (last level L) grid if(\(l\)+\(1==L\)) then // Direct solver on CPU \(\mathbf{s}_{1+1}\leftarrow\mathbf{K}_{1+1}\)-\(\mathbf{r}_{1+1}\) ; else \(\mathbf{s}_{1+1}\leftarrow\) V-cycle(\(\mathbf{r}_{1+1}\), 0,1+1) ; // Recursion \(\mathbf{s}_{1}\leftarrow\mathbf{s}_{1}+P_{1}\) ; // Prolongation // Post-smooth/-relaxation: smooth \(\mu_{2}\) times on \(\mathbf{K}_{1}\)\(\mathbf{s}_{1}=\mathbf{n}_{1}\) \(\mathbf{s}_{1}\leftarrow\)\(\mathbf{Smooth}(\mathbf{r}_{1}\),\(\mathbf{K}_{1}\),\(\mathbf{s}_{1}\),\(\mu_{2})\) ;  end if
```  
**Algorithm 2** AMG preconditioner (V-cycle). ```
Input:\(\mathbf{n}_{1},\mathbf{s}_{1},1\) Output:\(\mathbf{s}_{1}\) // Pre-smooth/-relaxation: smooth\(\mu_{1}\) times on \(\mathbf{K}_{1}\)\(\mathbf{s}_{1}=\mathbf{r}_{1}\) \(\mathbf{s}_{1}\leftarrow\)\(\mathbf{Smooth}(\mathbf{r}_{1}\),\(\mathbf{K}_{1}\),\(\mathbf{s}_{1}\),\(\mu_{1})\) ; \(\mathbf{r}\leftarrow\mathbf{r}\) - \(\mathbf{K}\) ; \(\mathbf{s}_{1}\leftarrow\)\(\mathbf{Smooth}(\mathbf{r}_{1}\),\(\mathbf{K}_{1}\),\(\mathbf{s}_{1}\),\(\mu_{2})\) ;  end if
```  
**Algorithm 3** Conjugate gradient algorithm  
We have several choices for the smoothing or relaxation operator mentioned in Algorithm 2. We can mention Gauss-Seidel, damped Jacobi, sparse approximate inverse (SPAI), and incomplete LU factorization multigrid smoothers, to name but a few. In the numerical experiments of this work, we adopt the simple and parallel smoother SPAI-0 [57] and the incomplete LU factorization with zero fill-ins (ILU0).