Parallel Energy-Minimization Prolongation for Algebraic Multigrid

With the increasing availability of powerful computational resources, scientific and engineering applications are becoming more demanding in terms of both memory and CPU time. For common methods used in the numerical approximation to partial differential equations (e.g., finite difference, finite volume, or finite element), the resulting approximation can easily grow to several millions or even billions of unknowns. The efficient solution to the associated sparse linear system of equations

either as a stand-alone system or as part of a nonlinear solve process, often represents a significant computational expense in the numerical application. Thus, research on sparse linear solvers continues to be a key topic for efficient simulation at large scales. One of the most popular sparse linear solvers is algebraic multigrid (AMG) [5, 6, 27] because of its potential for $O(n)$ computational cost in number of degrees-of-freedom $n$ for many problem types.

A fast converging AMG method relies on the complementary action of relaxation (e.g., with weighted Jacobi) and coarse grid correction, which is a projection step focused on eliminating the error that is not reduced by relaxation. Even in a purely algebraic setting, the main algorithmic decisions in multigrid are often based on heuristics for elliptic problems. As a result, for more complex applications, traditional methods often break down, requiring additional techniques to improve accuracy with a careful eye on overall computational complexity.

Even with advanced AMG methods, robustness remains an open problem for a variety of applications, especially in parallel. Yet, there have been several advances in recent years that have significantly improved convergence in a range of settings. Adaptive AMG [19] and adaptive smoothed aggregation [9] are among early attempts to assess the quality of the AMG setup phase during the setup process, with the ability to adaptively improve the interpolation operators. Later works focus on extending the adaptive ideas to more general settings [20], and in particular, Bootstrap AMG [3] further develops the idea of adaptive interpolation with least-squares interpolation coupled with locally relaxed vectors and multilevel eigenmodes. Other advanced approaches have a focus on specific AMG components, such as energy minimization of the interpolation operator [21, 33, 28, 26, 23], generalizing the strength of connection procedure [25, 4], or by considering the nonsymmetric nature of the problem directly [24, 22].

While AMG robustness and overall convergence has improved with combinations of the advances above, the overarching challenge of controlling cost is persistent. In this paper, we make a number of related contributions with a focus on AMG effectiveness and efficiency at large scale. Our key contributions are as follows:

• •

  The quality and sparsity of tentative interpolation is improved through a novel utilization of sparse $QR$ and a new process for sparsity pattern expansion that targets locally full-rank matrices for improved mode interpolation constraints;

• •

  We accompany the energy minimization construction of interpolation with new energy and convergence monitoring, thus limiting the total cost;

• •

  We apply a new preconditioning technique for the energy minimization process based on Gauss-Seidel applied to the blocks;

• •

  We present the non-trivial and efficient parallel implementation in detail; and

• •

  We demonstrate improved convergence and computational complexity with several large scale experiments.

The remainder of the paper is as follows. We begin with the basics of AMG in Section 2. In Section 3, we derive the energy minimization process based on QR factorizations and introduce a method for monitoring reduction of energy in practice. Finally, we conclude with several numerical experiments in Section 5 along with a discussion on performance.

The effectiveness of AMG as a solver depends on the complementary relationship between relaxation and coarse-grid correction, where the error not reduced by relaxation on the fine grid (e.g., with weighted-Jacobi or Gauss-Seidel) is accurately represented on the coarse grid, where a complementary error correction is computed. For a more in-depth introduction to AMG, see the works [10, 29]. Here, we focus our description of AMG on the coarse grid and interpolation setup, which are most relevant to the rest of the paper.

Constructing the AMG coarse grid begins with a partition of the $n$ unknowns of $A$ into a C-F partition of $n_{f}$ fine nodes and $n_{c}$ coarse nodes: $\{0,\ldots,n-1\}=\mathcal{C}\cup\mathcal{F}$. From this, we assume an ordering of $A$ by F-points followed by C-points:

where for example, $A_{ff}$ corresponds to entries in $A$ between two F-points. We also assume $A$ is SPD so that $A_{cf}=A_{fc}^{T}$. In classical AMG, prolongation takes the form

where $W$ must be sparse (for efficiency) and represents interpolation from the coarse grid to fine grid $F$-points.

In constructing prolongation of the form Eq. 3, there are two widely accepted guidelines, the so-called ideal [8, 35] and optimal [34, 7] forms of prolongation. Although both of these are not feasible in practical applications, leading to very expensive and dense prolongation operators, the concepts behind their definition are valuable guides for constructing effective $P$.

Ideal prolongation is constructed by starting with the above C-F partition and constructing $P_{\text{id}}$ as

Making $P_{\text{id}}$ the goal for interpolation is motivated by Corollary 3.4 from the theoretical work [12]. Here, the main assumption is a classical AMG framework where $P$ is of the form in equation (3).¹¹1The other assumptions are specific choices for the map to F-points $S=[I,0]$, for the map to C-points $R=[0,I]$, and for relaxation $X=\|A\|I$. In this setting, the choice of $W=-A_{ff}^{-1}A_{fc}$ minimizes the two-grid convergence of AMG relative to the choice of $P$, i.e., relaxation is fixed. Motivating our later energy minimization approach, $P_{\text{id}}$ can be viewed as having zero energy rows, as $AP_{\text{id}}$ is zero at all F-rows. Additionally, this classical AMG perspective likely makes the task of energy minimization easier, in that the conditioning of $A_{FF}$ is usually superior to that of $A$.

With optimal interpolation, the goal of interpolation is to capture the algebraically smoothest modes in span$(P)$, i.e., the modes left behind by relaxation. More specifically following [7], let $\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{n}$ and $v_{1},v_{2},\dots,v_{n}$ be the ordered eigenvalues and eigenvectors of the generalized eigenvalue problem $Ax=\lambda\tilde{M}x$, where $\tilde{M}$ is the symmetrized relaxation matrix, e.g., the diagonal of $A$ for Jacobi. (See the work [7] for more details.) Then, the two-grid convergence of AMG is minimized if

Note, that no assumptions on the structure of $P$ are made, as in equation (3). Motivating our later energy minimization approach, equation (5) indicates that span$(P)$ should capture low-energy modes relative to relaxation, which our Jacobi or Gauss-Seidel preconditioned energy minimization approach will explicitly target. Moreover, our energy minimization approach will incorporate constraints which explicitly force certain vectors to be in span$(P)$, where these vectors are chosen to represent the $v_{i}$ with smallest eigenvalues.

The idea of energy minimization AMG with constraints has been exploited for both symmetric and non-symmetric operators in several works [21, 33, 28, 26, 23], and, though requiring more computational effort than classical interpolation formulas, often provides improved preconditioners that balance the extra cost.

The energy minimization process combines the key aspects of ideal and optimal prolongation. To define this, we first introduce $V$, a basis for the near kernel of $A$ or the lowest energy modes of $A$. Then, energy minimization seeks to satisfy two requirements:

1. 1.

   Range: The range of prolongation must include the near kernel $V$:

   (6)

   $$V\subseteq\mathrm{range}(P)$$

2. 2.

   Minimal: The energy of each column of $P$ is minimized:

   (7)

   $$P=\operatorname*{argmin}_{P}\left(\mbox{tr}(P^{T}AP)\right).$$

To construct a $P$ that contains $V$ in the range and has minimal energy, we next introduce the key components needed by most energy minimization approaches, namely

1. i)

   a sparsity pattern for efficient application of $P$;

2. ii)

   a constraint to enforce Eq. 6; and

3. iii)

   an approximate block diagonal linear system for solving equation Eq. 7.

For practical use in AMG, the prolongation operator must be sparse, therefore construction begins by defining a sparse non-zero pattern for $P$. Assume that a strength of connection (SoC) matrix $S$ is provided where nonzero entries denote a strong coupling between degrees-of-freedom [27, 25, 4]. Next, let $P_{0}$ be a tentative prolongation with non-zero pattern $\mathcal{P}_{0}$, to be used as an initial guess for $P$.²²2See Section 3.2 for our contributions regarding the construction of $P_{0}$, which is based on the adaptive algorithm [15]. $P_{0}$ can be defined similarly to the tentative prolongation from smoothed aggregation AMG [31] in that $P_{0}$ interpolates the basis $V$, but needs further improvement, for example, with energy minimization. We next obtain an enlarged sparsity pattern $\mathcal{P}$ by growing $\mathcal{P}_{0}$ to include all strongly connected neighbors up to distance $k$. Denoting with $\overline{P}$ the unitary matrix obtained from $P$ by replacing its non-zeros with unitary entries, this is equivalent to

where $\mathcal{P}$ is the pattern of $\overline{P}_{k}$ (see [26]).

For a constraint condition satisfying Eq. 6, we start by splitting the near kernel basis $V$ with the same C-F splitting as in Eq. 2:

Recalling the form of interpolation Eq. 3, the near kernel requirement for $P$ becomes

which is a set of $n_{f}$ conditions on the rows of $W$. By denoting $\noindent w_{i}^{T}$ as the $i$-th row of $W$ (and $P$) and $\noindent v_{i}^{T}$ as the $i$-th row $V$, condition Eq. 10 is then exploited row-wise as

Using the sparsity pattern $\mathcal{P}$, we rewrite Eq. 11 for only the nonzeros in each row $\noindent w_{i}$. Letting the index set $\mathcal{J}_{i}$ be the column indices in the $i$-th row of $P$, this becomes

where $\overline{\noindent w_{i}}=\noindent w_{i}(\mathcal{J}_{i})$ collects only the nonzeros of $\noindent w_{i}$. It is important to note that for each of the $n_{f}$ fine points, the constraints Eq. 12, are independent of each other, because each entry of $W$ appears in only one constraint. Denote by $\widetilde{p}$ the column vector collecting the nonzero entries of $P$ row-wise. Similarly, we denote by $p$ the column vector containing the nonzero entries of $P$ column-wise. By definition, $\widetilde{p}$ and $p$ have the same size, equal to the number of nonzeros in $P$. Next, we write the constraints in the following matrix form:

where $\widetilde{g}$ collects the $\noindent v_{i}$ in Eq. 12 into a single vector, and where $\widetilde{B}^{T}$, is a block diagonal matrix composed of $V_{c}(\mathcal{J}_{i},:)^{T}$ due to the independence of constraints.

Lastly, we describe the reduced linear system framework for approximating Eq. 7. Minimizing Eq. 7 is equivalent to minimizing the energy of the individual columns of $P$ on the prescribed nonzero pattern:

where $p_{i}$ is the $i$-th column of $P$ and $\mathcal{P}_{i}$ is the sparsity pattern of column $i$.

Next, let $\mathcal{I}_{i}$ be the set of nonzero row indices of the $i$-th column of $W$ and let $\overline{h}_{i}$ be the vector collecting the nonzero entries of the $i$-th column of $W$. Then the minimization in Eq. 14 defines $\overline{h}_{i}$ with

where $A(\mathcal{I}_{i},\mathcal{I}_{i})$ is a square, relatively dense, submatrix of $A$ corresponding to the allowed nonzero indices $\mathcal{I}_{i}$ and $A(\mathcal{I}_{i},i)$ is a vector corresponding to the $i$-th column of $A$ at the allowed nonzero indices. Also in this case, each column of $P$ satisfies Eq. 14 independently. Thus, denoting by $p$ the column vector collecting the nonzero entries of $P$ column-wise (i.e., a rearranged version of $\widetilde{p}$) and denoting by $f$ the vector collecting each $-A(\mathcal{I}_{i},i)$, the minimization Eq. 14 is recast as

where $K$ is block diagonal with $A(\mathcal{I}_{i},\mathcal{I}_{i})$ on the $i$-th block.

The two conditions, one on the range of $P$ and one on the minimality of $P$, together form a constrained minimization problem, whose solution is the desired energy minimal prolongation. Casting this problem using Lagrange multipliers results in the saddle point system

The elements in Eq. 17 are the same as those defined in Eqs. 13 and 16, with the exception that $\widetilde{B}$, $\widetilde{g}$, and $\widetilde{p}$ are reordered following the columns of $P$, and $\lambda$ is the vector of Lagrange multipliers whose values are not needed for the purpose of setting up the prolongation. We emphasize that, with the entries of $p$ enumerated column-wise with respect to $P$, $K$ is block diagonal. Likewise, if $p$ is enumerated following the rows of $P$, then $B$ becomes block diagonal. Unfortunately, there is no sorting of $p$ able to make both $K$ and $B$ block diagonal at the same time. Nevertheless, it is possible to take advantage of this underlying structure in the numerical implementation, as will be shown later. Leveraging the block structure of $B$ is also important, because, as we will see in Section 3.1, our algorithm to minimize energy requires several applications of the orthogonal projector given as

The system Eq. 17 follows closely the method from  [26]. In sections 3.2–3.4, we outline our proposed improvements to energy minimization.

Energy minimization has historically been considered computationally expensive for AMG setup and real applications. Nevertheless, a cost effective implementation is still possible, but requires special care with the algorithm parallelization. In this work, we build our AMG preconditioner with Chronos [15, 13], which provides the standard numerical kernels usually required in a distributed memory AMG implementation such as the communication of ghost unknowns, coarse grid construction (e.g., computing the PMIS C-F partition [11]), sparse matrix-vector, and sparse matrix-matrix product, etc. For energy minimization, however, we developed three specific kernels that are not required in other AMG approaches, but are critical for an efficient parallel implementation here, i.e., the sparse product between $K$ and $p$, application of the projection $\Pi_{B}$, and symmetric Gauss-Seidel with $K$. We do not list Jacobi preconditioning, because it simply consists of a row-wise scaling of $P$.

The first issue related to the product by $K$ is that in practical applications $K$ cannot be stored. In fact, if we consider a prolongation $P$ having $r$ non-zeroes per row on average, the number of non-zeroes per column will be approximately $s\simeq\frac{n_{f}}{n_{c}}r$ and the $K$ matrix would be of size $n_{c}s\times n_{c}s$ with about $n_{c}s\,t\simeq n_{f}r\,t$ non-zeroes to be stored, where $t$ is the average number of non-zeroes per row of $A$. Often, $r\simeq t$ and storing $K$ becomes several times more expensive than storing $A$. For instance in practical elasticity problems, $K$ can be up to 20 times larger than $A$, making it unavoidable to proceed in matrix-free mode for $K$. Fortunately, the special structure of $K$ allows us to interpret the product $Kp$ as a sparse matrix-matrix (SpMM) product between $A$ and $P$, but with a fixed, prescribed pattern on $P$. This property can be easily derived from the definition of $K$ Eq. 16 and the vector form of the prolongation $p$. One advantage of prescribing a fixed pattern for $P$ is that the amount of data exchanged in SpMM is greatly reduced. First of all, the sparsity pattern adjacency information of $P$ can be communicated only once before entering PCG, then for all the successive $AP$ products only the entries of $P$ are exchanged. Moreover, all the buffers to receive and send messages through the network, can be allocated and set-up only once at the beginning and removed at the end of the minimization process. In practice, we find that these optimizations reduce the cost of SpMM by about 50%.

By contrast, the construction and application of $\Pi_{B}$ does not require any communication. In fact, we distribute the prolongation row-wise among processes, and since $B$ is block diagonal, with each block corresponding to a row of $P$ as shown in (13), we distribute $B$ blocks to the process owning the respective row of $P$. Following this scheme, each block $\mathbb{B}$ of $B$ is factorized independently to obtain $Q$ which is efficiently processed by the LAPACK routine DGEMV when applying $\Pi_{B}$.

Finally, we emphasize that there is no parallel bottleneck when parallelizing the application of symmetric Gauss-Seidel with $K$. This is because $K$ is block diagonal and, at least in principle, each block can be assigned to a different process. However as in the $Kp$ product, the main issue is the large size of $K$ which prevents explicit storage. As above, we rely on the equivalence between the $Kp$ product and the $AP$ product with a prescribed sparsity pattern. The symmetric Gauss-Seidel step is then performed matrix-free, similar to the $AP$ product, by saving again a considerable amount of data exchange because it is not necessary to communicate adjacency information or reallocate communication buffers. As expected, the symmetric Gauss-Seidel application exhibits a computational cost comparable to that of the $AP$ product.

In this section, we investigate the improved convergence offered by energy minimization and the resulting computational performance and parallel efficiency for real-world applications. This numerical section is divided into three parts: a detailed analysis on the prolongation energy reduction and related AMG convergence for two medium-size matrices, a weak scalability test for an elasticity problem on a uniformly refined cube, and a comparison study for a set of large real world matrices representing fluid dynamics and mechanical problems.

As a reference for the proposed approach, we compare with the well-known and open source solver GAMG [1], a smoothed aggregation-based AMG code from PETSc. In both cases, Chronos and GAMG are used with preconditioned conjugate gradients (PCG). The GAMG set-up is tuned for each problem starting from its default parameters, as suggested by the user guide [30]. As smoother we have chosen the most powerfull option available in Chronos and GAMG, that is FSAI and Chebyshev accelerated Jacobi, respectively. Each time such default parameters are modified, we report the chosen values in the relevant section. Since all of the experiments are run in parallel, the system matrices are partitioned with ParMETIS [18] before the AMG preconditioner set-up to reduce communication overhead.

All numerical experiments have been run on the Marconi100 supercomputer located in the Italian consortium for supercomputing (CINECA). Marconi100 consists of 980 nodes based on the IBM Power9 architecture, each equipped with two 16-core IBM POWER9 AC922 @3.1 GHz processors. For each test, the number of nodes, $N$, is selected so that each node has approximately 640,000,000 nonzero entries, and, consequently, the number of nodes is problem dependent. Each node is always fully exploited by using 32 MPI tasks, i.e., each task (core) has an average load of 20,000,000 nonzero entries. The number of cores is denoted $n_{cr}$. Only during the smaller cases in the weak scalability analysis are nodes partially used (i.e., with less than 32 MPI tasks). Even though Chronos can exploit hybrid MPI-OpenMP parallelism, for the sake of comparison, it is convenient to use just one thread, i.e., pure MPI parallelism. Moreover, for such a high load per core, we do not find that fine-grained OpenMP parallelism is of much help.

The numerical results are presented in terms of total number of computational cores used, $n_{cr}$, the grid and operator complexities, $C_{gd}$ and $C_{op}$, respectively, the number of iterations, $n_{it}$, and the setup, iteration, and total times, $T_{p}$, $T_{s}$, and $T_{t}$ = $T_{p}$ + $T_{s}$, respectively. For all the test cases, the right-hand side is a unit vector. The linear systems are solved with PCG and a zero initial guess, and convergence is achieved when the $\ell_{2}$-norm of the iterative residual drops below 8 orders of magnitude with respect to the $\ell_{2}$-norm of the right-hand side.