2.1 Test Problem: Poisson Equation

The Poisson problem was selected to test the performance of the implemented multigrid method. Poisson equation is an elliptic partial differential equation with several applications in science, which is commonly used for the evaluation of multigrid solvers (see the review of related works in Section 3).

The considered Poisson problem is defined on the unit square as domain, and Dirichlet boundary conditions are assumed, as expressed in Eq. 1. There, Δu(x,y)=uxx(x,y)+uyy(x,y) is the Laplacian operator; v(x,y) is a known source function; and g(x,y) is also a known function that determines the boundary conditions. The Poisson equation is discretized by using centered finite differences, in a regular grid of n+1 vertices per dimension. The discretization generates a sparse linear system, with (n−1)2 unknowns (one for each interior vertex) [²]. The resulting linear system is solved with multigrid to compute an estimation of u(x,y) in the points of the grid:

2.2 Multigrid Methods

Multigrid methods work under the idea of accelerating an iterative solver by using information provided by a global correction procedure, which operates on coarse grids to solve a fine grid problem which in turn is easier to solve. The recursive process, called the multigrid cycle, is applied until a direct solver can be applied without additional computation cost on the coarsest grid. V-cycle is one of the most popular type of multigrid cycles, characterized by computing the coarse grid correction down levels (the restriction phase), until finding a level where the direct method is applied.

Numerical computations on steps performed in coarse grids are quicker, as few vertices are considered, and the numerical error converges faster as the method goes down levels [⁸]. Then, an interpolation is applied (the prolongation phase) to determine values on upper level grids, until the finer grid is again reached. Multigrid methods are recursive in nature and a V-cycle can be extended to any number of levels.

In this article, geometric multigrid, implemented with cycles of type V, is used as a solver for the discretized Poisson problem. The problem domain is first partitioned into a uniform grid, in which an initial estimation is find by applying few iterations from an iterative smoother method. The residual of this initial estimation is then restricted into a a coarser grid, where it is corrected by applying few iterations of the smoother.

This process is repeated by taking successively coarser grids. When the coarsest grid is reached, estimations from the different grid levels are combined by interpolation to the finer grids, until the initial grid is reached. This process is illustrated in Fig. 1, where h is the width of the initial grid, n is the number of unknowns in each dimension of this grid, and R and I represent the restriction and interpolation operations, respectively.

Fig. 1 A non-truncated V cycle for MG 

A recursive implementation for a multigrid V cycle is presented in Algorithm 1, following the idea by Briggs et al. [²]. The implementation uses Gauss-Seidel (GS) as smoother method. GS applies a Red-Black update order strategy and executes a fixed number of iterations θ. Rh2h represents the restriction of the grid vertices to the next coarser grid; which is done by full-weighting of adjacent vertices. I2hhdenotes the interpolation to the finer grid, which is performed by bi-linear interpolation of adjacent vertices. Finally, the Successive Over-Relaxation (SOR) method is applied as the coarse grid solver.

Algorithm 1 Recursive multigrid V cycle: Vh(Ah,rh,u0) 

Gauss-Seidel Red-Black (GS-RB) is selected as smoother method, since it is a more parallelizable version of the traditional Gauss-Seidel solver. GS-RB is obtained by modifying the order in which grid vertices are updated. Vertices are first separated into two colors: red and black, intercalated. Then, they are updated with the usual GS expression, first applied to all red vertices, and then to the black ones [¹⁵]. This way, vertices of the same color may be updated concurrently. An important part of the smoother is the residual, as it is used as the independent term of the linear system solved by the multigrid method in each grid. In the case of GS-RB, the residual for black vertices is always null, as they are updated last. For the red vertices update, the residual coordinates associated to uk+1, are given by Eq. 2 (for i+j=even):

When solving a discretized Poisson linear system, GS-RB requires O(n4) operations in order to estimate a solution satisfying ‖rk‖<ϵ‖r0‖, ϵ>0. More efficient methods are Conjugate Gradient or Successive Over-Relaxation (SOR); both with O(n3) operations [⁴]. For the proposed implementation, SOR was chosen as coarse grid solver, as it may be obtained easily from the GS implementation. Specifically, SOR with Red-Black order (SOR-RB) is used. The linear system of interest is symmetric and positive definite, which implies that SOR-RB is convergent for any w∈(0,2) [⁴]. The value of w was chosen to maximize the convergence speed, as wopt=2/(1+sin⁡ (πh)). The coordinates of SOR residual, associated to uk+1, are given by Eq. 3 (valid for red and black vertices):

Each multigrid iteration applies one V cycle, with initial condition taken as the previous V cycle final estimation: uk+1=Vh(A,b,uk).

When developing a parallel implementation of a multigrid method in a distributed memory system, the number of vertices assigned to each processing unit decreases exponentially while descending in a V cycle. Thus, the cost of message passing may start to prevail over the cost of computations, producing a degradation of the overall computational efficiency. A technique to overcome this problem is, instead of reaching the coarsest grid, using truncated V cycles [¹⁶], as illustrated in Fig. 2. When using a truncated V cycle, Algorithm 1 iterates until reaching the coarsest grid, which is determined by the truncation criteria.

Fig. 2 A truncated V cycle for MG 

The truncated V cycle strategy has a specific drawback: if the cycle is truncated too soon, the (truncated) coarse grid may have too many vertices, and solving the coarse grid linear system may become a bottleneck, affecting the overall computational efficiency of the solver. Thus, selecting an appropriate truncation level that takes into account the resulting trade-off is one of the most important challenges when implementing distributed memory multigrid methods [¹⁰].

4.1 Truncated V Cycle: Criteria and Properties

The proposed method is based on choosing a truncation level which is deep enough to have few vertices per processing unit, but not as few as to have idle processing units. The formal definition is presented next.

Definition 1 (Truncation criteria). Consider an initial grid with n=q×2r vertices per dimension, where q is odd. In this case, each V cycle has a maximum of r levels. Given N≥1 and p processing units, truncate the V cycle at the deepest possible level k≤r, which also satisfies: nl2/p≥N, where nl=n/2l is the number of vertices per dimension in the (truncated) coarse grid.

Algorithm 2 presents a pseudocode for the proposed truncation criteria. Two conditions (line 3) control the iterative descending procedure. When using this criteria with an ideal load balancing, no processing unit is idle at the coarsest (truncated) grid. To compensate for some unbalance, the value of N should be chosen greater than one.

Algorithm 2 Proposed criteria for truncating a V cycle of geometric multigrid 

The level reached by the proposed criteria depends on the initial size n, the number of processing units p, and the value chosen for N. Theorem 1 characterizes this dependence.

Theorem 1. Let assume that the initial number of vertices per dimension is n=q×2r, q odd, r∈ℕ, r≥0. Then, the level l of the truncated V cycle reached with the proposed criteria, is (for p∈ℕ, p≥1):

Proof.

Example 1. Consider n=10240=5×211=q×2r and p=64=26. Without truncating the V cycle, the deepest possible level is l=11; resulting in a coarse grid with nc=5 vertices per dimension, and nc2/p=25/64<1 vertices per processing unit (for an ideal load balance). Now consider truncating each cycle with the proposed criteria, using N=4. In this case: n/Np=5×27<211. Therefore, by Theorem 1, the value of l is given by Eq. 4:

The number of vertices per dimension at this level is nl=n/29=20. Assuming an ideal load balance, the number of vertices per processing unit is nl2/p=6.25. Fig. 3 shows the reached level l, coarse grid vertices per dimension nl, and per processing unit nl2/p, for different values of p and N, as given by Theorem 1. Corollary 1 bounds the number of vertices in the coarse grid, independent from n and p.

Fig. 3 Truncated V cycle for n=5×211, when using the proposed criteria: level reached (l), coarse grid vertices per dimension (nl), and per processing unit nl2/p 

Corollary 1. For a truncated coarse grid reached at level l∈{1,…,r−1}, the total number of vertices is bounded by Np≤nc2<4Np, and the number of vertices per processing unit is bounded by N≤nc2/p<4N (for ideal load balance).

Proof. From the proof of Theorem 1:

n/2Np<2l≤n/Np⇔Np≤nc=n/2l<2Np⇔Np≤nc2<4Np.

4.2 Design and Implementation Details

This subsection presents the design considerations and implementation details of the proposed method.

4.2.2 Domain Decomposition and Processing

The decomposition approach is based on dividing the domain grid into rectangular sub-domains, to be assigned to each of the available processing units. A Cartesian grid is created and the optimal number of processing units per dimension (Nx and Ny) is determined using the MPI_Dims_create function of MPI. The output information is then used to create a new MPI communicator with Cartesian topology, by using the MPI_Cart_create function.

The grid vertices are then distributed along the processing units. Let (qx,qy) the coordinates of a generic processing unit q in the Cartesian communicator, processing unit q receives the inner grid vertices with coordinates in the two intervals Ixq and Iyq, as defined in Eqs. 5 and 6, where the integer division operator is used:

Ixq=[qx(n−1)Nx]+1, (qx+1)(n−1)Nx, (5)

Iyq=[qy(n−1)Ny]+1, (qy+1)(n−1)Ny. (6)

Once the data is partitioned, each processing unit applies the same MPI multigrid code to its own assigned vertices. The numerical resolution is performed in a coordinate way with other processing units, by exchanging halo layers of ghost vertices, as explained in the next subsection. Halo layers act as boundary values for each sub-domain, and contain values estimated by adjacent processing units. Fig. 4 illustrates an example of the domain decomposition.

Fig. 4 Example of domain decomposition into four sub-domains, with an overlap region (halo layers) of width one. Inner vertices (white circles), boundary vertices (white squares), and ghost vertices (black circles/squares) are distinguished 

5.1 Development and Execution Platform

The proposed multigrid method was implemented in the C programming language, using the MPICH implementation of the MPI standard (http://www.mpich.org). All floating point calculations use double precision.

The experimental evaluation was performed on HP ProLiant DL380 G9 servers with two Xeon Gold 6138 processors with 20 cores at 2.00 GHz each, and 128 GB of RAM, connected by Ethernet 10 Gbps, and running the Linux CentOS 7 OS, from National Supercomputing Center (Cluster-UY), Uruguay [¹³].

5.2 Test Problem, Convergence Criterion and Parameters

Tests were performed considering the Poisson problem defined in Eq. 1, whose exact solution given by Eq. 7:

The source function is taken as v(x,y)=−Δu(x,y), and boundary values g(x,y) are obtained by restricting the exact solution to the domain boundary.

Regarding the convergence criterion, SOR-RB stops execution when the residual of the corresponding linear system at step k satisfies ‖rk‖​2/‖r0‖​2≤10−6. The same criteria is used for the multigrid main iteration, but with b instead of r0. For multigrid the maximum number of iterations is set to 15 V cycles. For SOR-RB this is set to 5nc iterations; where nc is the number of vertices per dimension in the coarse grid. The truncation criteria uses N=5.

The first multigrid iteration (first V cycle) uses u0=0→. For subsequent cycles, the initial estimation is taken as the final estimation of the previous cycle (“warm restart”). The GS-RB smoother does θ=2 iterations before descending to a coarser grid (pre-smooth), and after ascending to a finer grid (post-smooth).

5.3 Execution Time, Speedup and Efficiency

The algorithmic speedup metric is applied to analyzed the performance of the proposed method. Five independent executions were performed to reduce bias due to the utilization of a non-dedicated hardware platform. The chosen values for n are of the form q×2r, with q∈{1,3,5}.

Table 1 reports the average execution times of the proposed method (Tp¯), with their corresponding coefficient of variation (V^(%)), using an increasing number of computing units (p, up to 70) and cluster nodes (c), for different problem sizes (n from 8192 to 20480).

The execution time results reported in Table 1 Results demonstrate that the proposed implementation was able to steadily reduce the execution times of the multigrid method, as p increases. In particular, for n=20480, the execution time reduced from 781 seconds for the sequential implementation to 14.0 seconds with p=70. Despite using a non dedicated cluster, a robust behavior was observed in the execution times, with coefficients of variation less than 5% for most cases.

Table 2 reports the corresponding algorithmic speedup values (Sp) for the execution times reported in Table 1. A graphical analysis of speedup values is presented in Fig. 5.

Fig. 5 Algorithmic speedup of the multigrid method with truncated V cycle (N=5) 

The efficiency results in Table 2 indicate that the best speedup (S70=55.6) was obtained for the two largest problems. Constant speedup intervals are observed at different values of p; particularly for [48,50] and [56,60], possibly explained by the transition from 2 to 3 computing nodes, and the corresponding increase in communication latencies. In general, speedup is always increasing for the selected range of processing units, which may indicate that communications do not prevail over computations, which is the main objective of the proposed truncation level criteria.

Table 3 reports the computational efficiency of the proposed method, defined by the normalized value of the speedup when using p computing resources (Equation 8):

The computational efficiency values ranged from 0.9 to 0.6. From p=16 to p=70, efficiency remained almost constant, in the order of 0.8 for the two largest problems, and 0.7 for the rest.

The best efficiency was 0.88, obtained for n=20480, when using p=24 in a two nodes environment.

A graphical analysis of efficiency values is presented in Fig. 6.

Fig. 6 Efficiency of the proposed multigrid algorithm with truncated V cycle (N=5) 

5.4 Scalability Analysis

A scalability analysis was performed to determine the capability of solving problems in rather similar execution times.

The setup for the analysis consisted in increasing the number of processing units p, together with the problem size n, whereas keeping constant the initial number of vertices per processing unit.

The value n2/p=40962 was chosen to have p=1 for the smallest problem size, to fit in 1GB of RAM. Values of n are of the form q×2r, q∈{1,3}.

Table 4 reports the average execution times and coefficients of variation of the proposed multigrid implementation. Results demonstrate a good scalability behavior of the proposed implementation.

The execution time only increased 53% when the problem size n increased in a factor of 8 (from 4096 to 32768). Also, the execution time remained constant for n in 4096 to 12288 and 20480 to 32768.