# Multi-GPU Acceleration for Finite Element Analysis in Structural Mechanics
## 3 Aggregation Algebraic Multigrid (AMG)  
Consider a linear system of equations,  
\[Ku=f, \tag{1}\]  
where \(K\in\mathbb{R}^{n\times n}\) is the coefficient matrix of the linear elasticity system, \(f\in\mathbb{R}^{n\times 1}\) is the right-hand side vector, \(u\in\mathbb{R}^{n\times 1}\) is the solution vector, and \(n\) is the number of unknowns. Multigrid methods use a two-grid scheme to address the problem. Let \(l\) be the grid level and \(n_{l}\) the number of unknowns in such a \(l\) level. Considering the initialization \(K_{0}=K\) and \(n_{0}=n\), we define the coarse grid coefficient matrices \(K_{l+1}\) recursively as follows:  
\[K_{l+1}=R_{l}K_{l}P_{l}, \tag{2}\]  
where \(P_{l}\in\mathbb{R}^{n_{l}\times n_{l+1}}\) is the interpolation or prolongation operator, \(R_{l}\in\mathbb{R}^{n_{l+1}\times n_{l}}\) (typically obtained as \(P_{l}^{T}\)) is the restriction operator, and \(K_{l}\in\mathbb{R}^{n_{l}\times n_{l}}\) and \(K_{l+1}\in\mathbb{R}^{n_{l+1}\times n_{l+1}}\) are the fine and coarse grid coefficient matrices, respectively. We calculate these transfer operators recursively until a grid level \(L\), where the number of unknowns \(n_{L}\) is small enough to solve it within a reasonable time, usually using some direct solver. Algorithm 1 outlines the AMG setup stage for aggregation-based AMG methods [53], where we calculate the grid hierarchy and the transfer operators using the algebraic proper  
Figure 1: Multi-GPU node architecture. ties of the coefficient matrix \(K\). This stage requires the coefficient matrix \(K\in\mathbb{R}^{n\times n}\) and the near null-space vector \(B\in\mathbb{R}^{n\times m}\), with \(m\) number of vectors representing low eigenmodes of the problem. The setup phase requires the following operations: strength, aggregate, tentative, prolongate, and Galerkin projection. The strength operator constructs a graph \(C_{l}\) of strong connections at the corresponding level \(l\)[36; 54]. We build this operator using a coarsening ratio by evaluating the relationship between unknowns with the diagonal and outside diagonal coefficients. Using a block matrix by vertex (grouping into blocks all unknowns of a grid point) instead of considering each unknown as a vertex of \(C_{l}\) often shows computational advantages. ```
Input:\(\mathbf{K}\), \(\mathbf{B}\), \(n\) Output:\(\mathbf{K_{1}},\ldots,\mathbf{K_{n_{l}}}\), \(\mathbf{P_{0}},\ldots,\mathbf{P_{n_{l}-1}}\), \(\mathbf{R_{0}},\ldots,\mathbf{R_{n_{l}-1}}\) \(\mathbf{K_{0}}\leftarrow\mathbf{K}\), \(\mathbf{B_{0}}\leftarrow\mathbf{B}\), \(n_{0}\gets n\), \(l\leftarrow\emptyset\), \(L\gets l\) ; while\(n_{l}\) is too big to solve \(\mathbf{C_{l}}\leftarrow\mathit{strength}(\mathbf{K_{1}})\) ; // Strength-of-connections \(\mathbf{Agg_{i}}\leftarrow\mathit{aggregate}(\mathbf{C_{l}})\) ; // Construct coarse aggregates \(\tilde{\mathbf{P}}_{\mathbf{l}}\), \(\mathbf{B_{i+1}}\leftarrow\mathit{prolongate}(\mathbf{Agg_{i}},\mathbf{B_{i}})\) ; // Tentative \(\mathbf{P_{l}}\leftarrow\tilde{\mathbf{S_{i}}}\tilde{\mathbf{P_{l}}}\) ; // Improve tentative interpolation \(\mathbf{R_{l}}\leftarrow\mathbf{P_{l}^{T}}\); \(\mathbf{K_{l+1}}\leftarrow\mathbf{R_{l}}\mathbf{K_{l}}\mathbf{P_{l}}\) ; // Galerkin projection \(l\gets l+1\), \(L\gets l\), \(n_{0}\gets n\);  end while
```  
**Algorithm 1**AMG setup  
We then use the edges of the graph \(C_{l}\) for the aggregation process. The problem consists of selecting a set of independent grid points as root nodes, which we group with the neighbors to form the aggregates. This grouping can be performed in parallel using the Parallel Maximal Independent Set (PMIS) algorithm [55; 56]. Then, the tentative prolongation matrix \(\tilde{P}_{l}\) is built as a simple grid transfer operator by a piecewise constant interpolation, where each row corresponds to a grid point of the coarse grid and the columns correspond to an aggregate of \(Agg_{l}\).