Search

This paper presents performance results comparing MPI-based implementations of the popular Conjugate Gradient (CG) method and several of its communication hiding (or 'pipelined') variants. Pipelined CG methods are designed to efficiently solve SPD linear systems on massively parallel distributed memory hardware, and typically display significantly improved strong scaling compared to classic CG. This increase in parallel performance is achieved by overlapping the global reduction phase (MPI_Iallreduce) required to compute the inner products in each iteration by (chiefly local) computational work such as the matrix-vector product as well as other global communication. This work includes a brief introduction to the deep pipelined CG method for readers that may be unfamiliar with the specifics of the method. A brief overview of implementation details provides the practical tools required for implementation of the algorithm. Subsequently, easily reproducible strong scaling results on the US Department of Energy (DoE) NERSC machine 'Cori' (Phase I - Haswell nodes) on up to 1024 nodes with 16 MPI ranks per node are presented using an implementation of p(l)-CG that is available in the open source PETSc library. Observations on the staggering and overlap of the asynchronous, non-blocking global communication phases with communication and computational kernels are drawn from the experiments., Comment: EuroMPI 2019, 10 - 13 September 2019, ETH Zurich, Switzerland, 11 pages, 4 figures, 1 table. arXiv admin note: substantial text overlap with arXiv:1902.03100

Pipelined Krylov subspace methods (also referred to as communication-hiding methods) have been proposed in the literature as a scalable alternative to classic Krylov subspace algorithms for iteratively computing the solution to a large linear system in parallel. For symmetric and positive definite system matrices the pipelined Conjugate Gradient method outperforms its classic Conjugate Gradient counterpart on large scale distributed memory hardware by overlapping global communication with essential computations like the matrix-vector product, thus hiding global communication. A well-known drawback of the pipelining technique is the (possibly significant) loss of numerical stability. In this work a numerically stable variant of the pipelined Conjugate Gradient algorithm is presented that avoids the propagation of local rounding errors in the finite precision recurrence relations that construct the Krylov subspace basis. The multi-term recurrence relation for the basis vector is replaced by two-term recurrences, improving stability without increasing the overall computational cost of the algorithm. The proposed modification ensures that the pipelined Conjugate Gradient method is able to attain a highly accurate solution independently of the pipeline length. Numerical experiments demonstrate a combination of excellent parallel performance and improved maximal attainable accuracy for the new pipelined Conjugate Gradient algorithm. This work thus resolves one of the major practical restrictions for the useability of pipelined Krylov subspace methods., Comment: 15 pages, 5 figures, 1 table, 2 algorithms

Pipelined Krylov subspace methods avoid communication latency by reducing the number of global synchronization bottlenecks and by hiding global communication behind useful computational work. In exact arithmetic pipelined Krylov subspace algorithms are equivalent to classic Krylov subspace methods and generate identical series of iterates. However, as a consequence of the reformulation of the algorithm to improve parallelism, pipelined methods may suffer from severely reduced attainable accuracy in a practical finite precision setting. This work presents a numerical stability analysis that describes and quantifies the impact of local rounding error propagation on the maximal attainable accuracy of the multi-term recurrences in the preconditioned pipelined BiCGStab method. Theoretical expressions for the gaps between the true and computed residual as well as other auxiliary variables used in the algorithm are derived, and the elementary dependencies between the gaps on the various recursively computed vector variables are analyzed. The norms of the corresponding propagation matrices and vectors provide insights in the possible amplification of local rounding errors throughout the algorithm. Stability of the pipelined BiCGStab method is compared numerically to that of pipelined CG on a symmetric benchmark problem. Furthermore, numerical evidence supporting the effectiveness of employing a residual replacement type strategy to improve the maximal attainable accuracy for the pipelined BiCGStab method is provided., Comment: 27 pages, 9 figures, 1 table

Krylov subspace methods are widely known as efficient algebraic methods for solving large scale linear systems. However, on massively parallel hardware the performance of these methods is typically limited by communication latency rather than floating point performance. With HPC hardware advancing towards the exascale regime the gap between computation and communication keeps steadily increasing, imposing the need for scalable alternatives to traditional Krylov subspace methods. One such approach are the so-called pipelined Krylov subspace methods, which reduce the number of global synchronization points and overlap global communication latency with local arithmetic operations, thus hiding the global reduction phases behind useful computations. To obtain this overlap the traditional Krylov subspace algorithm is reformulated by introducing a number of auxiliary vector quantities, which are computed using additional recurrence relations. Although pipelined Krylov subspace methods are equivalent to traditional Krylov subspace methods in exact arithmetic, local rounding errors induced by the multi-term recurrence relations in finite precision may in practice affect convergence significantly. This numerical stability study aims to characterize the effect of local rounding errors on attainable accuracy in various pipelined versions of the popular Conjugate Gradient method. Expressions for the gaps between the true and recursively computed variables that are used to update the search directions in the different CG variants are derived. Furthermore, it is shown how these results can be used to analyze and correct the effect of local rounding error propagation on the maximal attainable accuracy of pipelined CG methods. The analysis in this work is supplemented by numerical experiments that demonstrate the numerical behavior of the pipelined CG methods., Comment: 28 pages, 7 figures, 1 table, 4 algorithms

Krylov subspace methods are among the most efficient solvers for large scale linear algebra problems. Nevertheless, classic Krylov subspace algorithms do not scale well on massively parallel hardware due to synchronization bottlenecks. Communication-hiding pipelined Krylov subspace methods offer increased parallel scalability by overlapping the time-consuming global communication phase with computations such as spmvs, hence reducing the impact of the global synchronization and avoiding processor idling. One of the first published methods in this class is the pipelined Conjugate Gradient method (p-CG). However, on large numbers of processors the communication phase may take much longer than the computation of a single spmv. This work extends the pipelined CG method to deeper pipelines, denoted as p(l)-CG, which allows further scaling when the global communication phase is the dominant time-consuming factor. By overlapping the global all-to-all reduction phase in each CG iteration with the next l spmvs (deep pipelining), the method hides communication latency behind additional computational work. The p(l)-CG algorithm is derived from similar principles as the existing p(l)-GMRES method and by exploiting operator symmetry. The p(l)-CG method is also compared to other Krylov subspace methods, including the closely related classic CG and D-Lanczos methods and the pipelined CG method by Ghysels et al.. By analyzing the maximal accuracy attainable by the p(l)-CG method it is shown that the pipelining technique induces a trade-off between performance and numerical stability. A preconditioned version of the algorithm is also proposed and storage requirements and performance estimates are discussed. Experimental results demonstrate the possible performance gains and the attainable accuracy of deeper pipelined CG for solving large scale symmetric linear systems., Comment: 34 pages (incl. 6 pages appendix), 12 figures, 2 tables, 5 algorithms

By reducing the number of global synchronization bottlenecks per iteration and hiding communication behind useful computational work, pipelined Krylov subspace methods achieve significantly improved parallel scalability on present-day HPC hardware. However, this typically comes at the cost of a reduced maximal attainable accuracy. This paper presents and compares several stabilized versions of the communication-hiding pipelined Conjugate Gradients method. The main novel contribution of this work is the reformulation of the multi-term recurrence pipelined CG algorithm by introducing shifts in the recursions for specific auxiliary variables. These shifts reduce the amplification of local rounding errors on the residual. The stability analysis presented in this work provides a rigorous method for selection of the optimal shift value in practice. It is shown that, given a proper choice for the shift parameter, the resulting shifted pipelined CG algorithm restores the attainable accuracy and displays nearly identical robustness to local rounding error propagation compared to classical CG. Numerical results on a variety of SPD benchmark problems compare different stabilization techniques for the pipelined CG algorithm, showing that the shifted pipelined CG algorithm is able to attain a high accuracy while displaying excellent parallel performance., Comment: 18 pages, 4 figures, 4 algorithms, 1 table

A High Performance Computing alternative to traditional Krylov subspace methods, pipelined Krylov subspace solvers offer better scalability in the strong scaling limit compared to standard Krylov subspace methods for large and sparse linear systems. The typical synchronization bottleneck is mitigated by overlapping time-consuming global communication phases with local computations in the algorithm. This paper describes a general framework for deriving the pipelined variant of any Krylov subspace algorithm. The proposed framework was implicitly used to derive the pipelined Conjugate Gradient (p-CG) method in "Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm" by P. Ghysels and W. Vanroose, Parallel Computing, 40(7):224--238, 2014. The pipelining framework is subsequently illustrated by formulating a pipelined version of the BiCGStab method for the solution of large unsymmetric linear systems on parallel hardware. A residual replacement strategy is proposed to account for the possible loss of attainable accuracy and robustness by the pipelined BiCGStab method. It is shown that the pipelined algorithm improves scalability on distributed memory machines, leading to significant speedups compared to standard preconditioned BiCGStab., Comment: 24 pages, 5 figures, 4 tables, 45 references

Pipelined Krylov subspace methods typically offer improved strong scaling on parallel HPC hardware compared to standard Krylov subspace methods for large and sparse linear systems. In pipelined methods the traditional synchronization bottleneck is mitigated by overlapping time-consuming global communications with useful computations. However, to achieve this communication hiding strategy, pipelined methods introduce additional recurrence relations for a number of auxiliary variables that are required to update the approximate solution. This paper aims at studying the influence of local rounding errors that are introduced by the additional recurrences in the pipelined Conjugate Gradient method. Specifically, we analyze the impact of local round-off effects on the attainable accuracy of the pipelined CG algorithm and compare to the traditional CG method. Furthermore, we estimate the gap between the true residual and the recursively computed residual used in the algorithm. Based on this estimate we suggest an automated residual replacement strategy to reduce the loss of attainable accuracy on the final iterative solution. The resulting pipelined CG method with residual replacement improves the maximal attainable accuracy of pipelined CG, while maintaining the efficient parallel performance of the pipelined method. This conclusion is substantiated by numerical results for a variety of benchmark problems., Comment: 26 pages, 6 figures, 2 tables, 4 algorithms

This paper introduces and explores the class of expansion preconditioners EX(m) that forms a direct generalization to the classic complex shifted Laplace (CSL) preconditioner for Helmholtz problems. The construction of the EX(m) preconditioner is based upon a truncated Taylor series expansion of the original Helmholtz operator inverse. The expansion preconditioner is shown to significantly improve Krylov solver convergence rates for the Helmholtz problem for growing values of the number of series terms m. However, the addition of multiple terms in the expansion also increases the computational cost of applying the preconditioner. A thorough cost-benefit analysis of the addition of extra terms in the EX(m) preconditioner proves that the CSL or EX(1) preconditioner is the practically most efficient member of the expansion preconditioner class. Additionally, possible extensions to the expansion preconditioner class that further increase preconditioner efficiency are suggested., Comment: 19 pages, 6 figures, 4 tables

This paper improves the convergence and robustness of a multigrid-based solver for the cross sections of the driven Schroedinger equation. Adding an Coupled Channel Correction Step (CCCS) after each multigrid (MG) V-cycle efficiently removes the errors that remain after the V-cycle sweep. The combined iterative solution scheme (MG-CCCS) is shown to feature significantly improved convergence rates over the classical MG method at energies where bound states dominate the solution, resulting in a fast and scalable solution method for the complex-valued Schroedinger break-up problem for any energy regime. The proposed solver displays optimal scaling; a solution is found in a time that is linear in the number of unknowns. The method is validated on a 2D Temkin-Poet model problem, and convergence results both as a solver and preconditioner are provided to support the O(N) scalability of the method. This paper extends the applicability of the complex contour approach for far field map computation [S. Cools, B. Reps, W. Vanroose, An Efficient Multigrid Calculation of the Far Field Map for Helmholtz and Schroedinger Equations, SIAM J. Sci. Comp. 36(3) B367--B395, 2014]., Comment: 24 pages, 10 figures, 1 table

Classical iterative methods for tomographic reconstruction include the class of Algebraic Reconstruction Techniques (ART). Convergence of these stationary linear iterative methods is however notably slow. In this paper we propose the use of Krylov solvers for tomographic linear inversion problems. These advanced iterative methods feature fast convergence at the expense of a higher computational cost per iteration, causing them to be generally uncompetitive without the inclusion of a suitable preconditioner. Combining elements from standard multigrid (MG) solvers and the theory of wavelets, a novel wavelet-based multi-level (WMG) preconditioner is introduced, which is shown to significantly speed-up Krylov convergence. The performance of the WMG-preconditioned Krylov method is analyzed through a spectral analysis, and the approach is compared to existing methods like the classical Simultaneous Iterative Reconstruction Technique (SIRT) and unpreconditioned Krylov methods on a 2D tomographic benchmark problem. Numerical experiments are promising, showing the method to be competitive with the classical Algebraic Reconstruction Techniques in terms of convergence speed and overall performance (CPU time) as well as precision of the reconstruction., Comment: Journal of Computational and Applied Mathematics (2014), 26 pages, 13 figures, 3 tables

In this paper we construct and analyse a level-dependent coarsegrid correction scheme for indefinite Helmholtz problems. This adapted multigrid method is capable of solving the Helmholtz equation on the finest grid using a series of multigrid cycles with a grid-dependent complex shift, leading to a stable correction scheme on all levels. It is rigourously shown that the adaptation of the complex shift throughout the multigrid cycle maintains the functionality of the two-grid correction scheme, as no smooth modes are amplified in or added to the error. In addition, a sufficiently smoothing relaxation scheme should be applied to ensure damping of the oscillatory error components. Numerical experiments on various benchmark problems show the method to be competitive with or even outperform the current state-of-the-art multigrid-preconditioned Krylov methods, like e.g. CSL-preconditioned GMRES or BiCGStab., Comment: 21 pages

In this paper we present a new highly efficient calculation method for the far field amplitude pattern that arises from scattering problems governed by the d-dimensional Helmholtz equation and, by extension, Schr\"odinger's equation. The new technique is based upon a reformulation of the classical real-valued Green's function integral for the far field amplitude to an equivalent integral over a complex domain. It is shown that the scattered wave, which is essential for the calculation of the far field integral, can be computed very efficiently along this complex contour (or manifold, in multiple dimensions). Using the iterative multigrid method as a solver for the discretized damped scattered wave system, the proposed approach results in a fast and scalable calculation method for the far field map. The complex contour method is successfully validated on Helmholtz and Schr\"odinger model problems in two and three spatial dimensions, and multigrid convergence results are provided to substantiate the wavenumber scalability and overall performance of the method., Comment: SIAM Journal on Scientific Computing, 29 pages, 10 figures, 5 tables

In this paper we solve the Helmholtz equation with multigrid preconditioned Krylov subspace methods. The class of Shifted Laplacian preconditioners are known to significantly speed-up Krylov convergence. However, these preconditioners have a parameter beta, a measure of the complex shift. Due to contradictory requirements for the multigrid and Krylov convergence, the choice of this shift parameter can be a bottleneck in applying the method. In this paper, we propose a wavenumber-dependent minimal complex shift parameter which is predicted by a rigorous k-grid Local Fourier Analysis (LFA) of the multigrid scheme. We claim that, given any (regionally constant) wavenumber, this minimal complex shift parameter provides the reader with a parameter choice that leads to efficient Krylov convergence. Numerical experiments in one and two spatial dimensions validate the theoretical results. It appears that the proposed complex shift is both the minimal requirement for a multigrid V-cycle to converge, as well as being near-optimal in terms of Krylov iteration count., Comment: 20 pages