We present a variant of the s-step Preconditioned Conjugate Gradient (PCG) method that combines a Chebyshev-stabilized Krylov basis with a Forward Gauss-Seidel (FGS) iteration for the solution of the reduced Gram systems. In s-step Conjugate Gradient, multiple search directions are generated per outer iteration, reducing global synchronization costs but requiring the solution of small dense Gram systems whose conditioning is critical for stability. We analyze the structure of the Chebyshev Gram matrix and show that its moment-based representation is associated with favorable conditioning properties for moderate step sizes. Building on inexact Krylov theory and on the classical equivalence between FGS and Modified Gram-Schmidt (MGS), we provide a structural analysis and theoretical rationale supporting the use of a small number of FGS sweeps, while preserving the convergence behavior observed in practical regimes. Large-scale experiments on modern NVIDIA GPU architectures demonstrate that the proposed Chebyshev-stabilized, Gauss-Seidel-enhanced s-step PCG achieves convergence comparable to classical CG while reducing synchronization overhead, making it a stable and scalable alternative for current and next-generation accelerator systems.

In this paper, we explore polynomial accelerators that are well-suited for parallel computations, specifically as smoothers in Algebraic MultiGrid (AMG) preconditioners for symmetric positive definite matrices. These accelerators address a minimax problem, initially formulated in Lottes (Numer. Lin. Alg. with Appl. 30(6), e2518 2023), aiming to achieve an optimal (or near-optimal) bound for a polynomial-dependent constant involved in the AMG V-cycle error bound, without requiring information about the matrices' spectra. Lottes focuses on Chebyshev polynomials of the 4th-kind and defines the relevant recurrence formulas applicable to a general convergent basic smoother. In this paper, we demonstrate the efficacy of these accelerations for large-scale applications on modern GPU-accelerated supercomputers. Furthermore, we formulate a variant of the aforementioned minimax problem, which naturally leads to solutions relying on Chebyshev polynomials of the 1st-kind as accelerators for a basic smoother. For all the polynomial accelerations, we describe efficient GPU kernels for their application and demonstrate their comparable effectiveness on standard benchmarks at very large scales.

The recycling of Krylov subspaces in iterative methods is a powerful strategy for reducing the computational burden of solving large-scale linear systems. However, traditional approaches often suffer from a loss of orthogonality among the basis vectors of the Krylov and recycled subspaces, leading to slower convergence rates. This paper introduces a novel hybrid Gauss--Seidel and Jacobi Gram--Schmidt algorithm that enhances both intra-and interblock orthogonalization, coupled with a block lower triangular correction matrix. Our approach ensures that Krylov vectors maintain orthogonality to working accuracy and demonstrates a significant reduction in the number of iterations required per linear system, outperforming established methods like HH-GMRES and the classical GCRO-DR. Additionally, we explore the utility of the recycle subspace eigenspectrum as an alternative to traditional preconditioners, with shifts that further decrease iteration counts. By employing the normwise relative backward error instead of the Arnoldi relative residual for convergence tests, our method avoids stagnation and enhances reliability. Applicable across a broad spectrum of recycling techniques, our algorithm offers substantial improvements in the efficiency and accuracy of solvers for large-scale linear systems.

Communication-avoiding Krylov methods require solving small dense Gram systems at each outer iteration. We present a low-synchronization approach based on Forward Gauss–Seidel (FGS), which exploits the structure of Gram matrices arising from Chebyshev polynomial bases. We show that a single FGS sweep is mathematically equivalent to Modified Gram–Schmidt (MGS) orthogonalization in the A-norm and provide corresponding backward error bounds. For weak scaling on AMD MI-series GPUs, we demonstrate that 20–30 FGS iterations preserve scalability up to 64 GPUs with problem sizes exceeding 700 million unknowns. We further extend this approach to Algebraic MultiGrid (AMG) coarse-grid solves, removing the need to assemble or factor dense coarse operators

We present SKA-SGD (Streaming Krylov-Accelerated Stochastic Gradient Descent), a novel optimization approach that accelerates convergence for ill-conditioned problems by projecting stochastic gradients onto a low-dimensional Krylov subspace. Directly inspired by recent advances in s-step Conjugate Gradient methods with streaming Gauss-Seidel Gram solvers \cite{dambra2025sstep}, our method extends these techniques to the stochastic optimization domain. Our approach combines three key innovations: (1) projection coefficients computed via a single streaming Gauss-Seidel iteration, which is mathematically equivalent to Modified Gram-Schmidt orthogonalization; (2) a Chebyshev polynomial basis for constructing the Krylov subspace, providing superior numerical stability; and (3) efficient implementation for AMD GPUs using HIP. We prove that our streaming approach achieves a backward error near machine precision with $O(s^2)$ complexity rather than $O(s^3)$, where $s$ is the Krylov subspace dimension. Experimental results demonstrate that SKA-SGD significantly outperforms standard SGD and Adam in convergence rate and final error, particularly for problems with condition numbers exceeding $10^3$. GPU performance analysis reveals a crossover point where communication-avoiding benefits outweigh computational overhead, typically occurring at moderate scale ($p \approx 64$ processors) for problem sizes $n \geq 10^6$.

Abstract. In this study, we introduce algorithms optimized for GPU architectures, aimed at efficiently solving large sparse linear systems, a central challenge in Navier–Stokes pressure projection problems. Our approach includes an adaptation of the GMRES algorithm, drawing inspiration from the merged vector operations first proposed by Bielich et al. [ Parallel Comput., 112 (2022), 102940]. This adaptation increases computational intensity on GPU platforms through optimized vector update strategies. The algorithm incorporates modified and classical Gram–Schmidt methods with an algebraic multigrid (AMG) preconditioner, each tailored for GPU performance. A key innovation in our work is the development of a Gram–Schmidt projector [Formula: see text] employing a rank-1 perturbation of the identity matrix. Designed to maximize the high memory bandwidth utilization of the AMD MI-250X GPU, this approach includes a strategy for treating the unit diagonal that minimizes memory reads, leading to a 25% increase in computational efficiency. The application of perturbation theory further ensures that orthogonality loss is limited to [Formula: see text], where [Formula: see text] is the number of iterations. Additionally, we introduce a mixed AMG [Formula: see text]-cycle strategy combining ILU(0) and [Formula: see text]-Jacobi smoothers, which achieves a 30–50% reduction in GPU compute times compared to conventional methods, while maintaining low backward error. This strategy, alongside our novel treatment of the diagonal in triangular matrices, marks a substantial increase in AMG efficicency for GPU systems. We believe that these contributions represent a significant advance in optimizing GMRES+AMG algorithms for GPU computations. The empirical results demonstrate notable speed increments and maintain rigorous backward error bounds, underscoring the potential of our methods to substantially increase computational efficiency in large-scale scientific applications.

Incomplete LU (ILU) smoothers are effective in the algebraic multigrid (AMG) V$$ V $$-cycle for reducing high-frequency components of the error. However, the requisite direct triangular solves are comparatively slow on GPUs. Previous work has demonstrated the advantages of Jacobi iteration as an alternative to direct solution of these systems. Depending on the threshold and fill-level parameters chosen, the factors can be highly nonnormal and Jacobi is unlikely to converge in a low number of iterations. We demonstrate that row scaling can reduce the departure from normality, allowing us to replace the inherently sequential solve with a rapidly converging Richardson iteration. There are several advantages beyond the lower compute time. Scaling is performed locally for a diagonal block of the global matrix because it is applied directly to the factor. Further, an ILUT Schur complement smoother maintains a constant GMRES iteration count as the number of MPI ranks increases, and thus parallel strong-scaling is improved. Our algorithms have been incorporated into hypre, and we demonstrate improved time to solution for linear systems arising in the Nalu-Wind and PeleLM pressure solvers. For large problem sizes, GMRES+$$ + $$AMG executes at least five times faster when using iterative triangular solves compared with direct solves on massively parallel GPUs. Convergence histories of the GMRES+AMG solver with a polynomial Guass-Seidel and the ILU direct and iterative smoothers, for a matrix from the projection step in PeleLM with 11M unknowns, solved on NREL's Eagle computer. The iterative ILU smoother shows a 5X speed-up over the direct smoother.image

The GMRES algorithm of Saad and Schultz (1986) is an iterative method for approximately solving linear systems $A{\bf x}={\bf b}$, with initial guess ${\bf x}_0$ and residual ${\bf r}_0 = {\bf b} - A{\bf x}_0$. The algorithm employs the Arnoldi process to generate the Krylov basis vectors (the columns of $V_k$). It is well known that this process can be viewed as a $QR$ factorization of the matrix $B_k = [\: {\bf r}_0, AV_k\:]$ at each iteration. Despite an ${O}(\epsilon)\kappa(B_k)$ loss of orthogonality, for unit roundoff $\epsilon$ and condition number $\kappa$, the modified Gram-Schmidt formulation was shown to be backward stable in the seminal paper by Paige et al. (2006). We present an iterated Gauss-Seidel formulation of the GMRES algorithm (IGS-GMRES) based on the ideas of Ruhe (1983) and \'{S}wirydowicz et al. (2020). IGS-GMRES maintains orthogonality to the level ${O}(\epsilon)\kappa(B_k)$ or ${O}(\epsilon)$, depending on the choice of one or two iterations; for two Gauss-Seidel iterations, the computed Krylov basis vectors remain orthogonal to working precision and the smallest singular value of $V_k$ remains close to one. The resulting GMRES method is thus backward stable. We show that IGS-GMRES can be implemented with only a single synchronization point per iteration, making it relevant to large-scale parallel computing environments. We also demonstrate that, unlike MGS-GMRES, in IGS-GMRES the relative Arnoldi residual corresponding to the computed approximate solution no longer stagnates above machine precision even for highly non-normal systems.

Optimizing a particular subprogram out of the set of Basic (sparse) Linear Algebra Subprograms (BLAS) for a given architecture is a common topic of research. In applications, however, these BLAS functions rarely appear in isolation; usually, many of them are used together, in various combinations and with varying inputs. As the need to solve a large, sparse linear system is ubiquitous throughout HPC applications, linear solvers constitute a realistic, sufficiently complex and well-defined representative use case for composite BLAS routines. To this end, based on a representative set of matrices drawn from a diverse set of fields, we present a framework to study, from the performance and energy perspective, the efficacy of GPU-resident parallel Conjugate Gradient (CG) linear solver with different preconditioner options, including Gauss-Seidel, Jacobi, and incomplete Cholesky. We also propose a novel GPU-based preconditioner, in which the triangular solves are approximated by an iterative process. The development of this preconditioner was motivated by solving large graph Laplacian linear systems, for which the existing preconditioners either perform slow on GPU-based platforms or are not applicable. We compare the performance of these preconditioners on different hardware accelerator architectures, i.e., AMD MI250X, MI100, Nvidia A100, V100, and Jetson. Our experiments reveal performance trade-offs and provide information on how to select the best strategy for the given linear system, dictated by its properties, and the platform of interest. We demonstrate the application of our novel preconditioner for solving CG and graph Laplacian systems. Overall, the framework can be utilized as a benchmark to guide informed decisions in choosing a specific preconditioner, i.e., whether it is better to rely on the performance of a triangular solver or on the performance of sparse matrix-vector product. Finally, by considering power consumption to solve the linear systems, we report the energy footprint for the solvers.

In this study, we introduce algorithms optimized for GPU architectures, aimed at efficiently solving large sparse linear systems, a central challenge in Navier--Stokes pressure projection problems. Our approach includes an adaptation of the GMRES algorithm, drawing inspiration from the merged vector operations first proposed by Bielich et al. [Parallel Comput., 112 (2022), 102940]. This adaptation increases computational intensity on GPU platforms through optimized vector update strategies. The algorithm incorporates modified and classical Gram--Schmidt methods with an algebraic multigrid (AMG) preconditioner, each tailored for GPU performance. A key innovation in our work is the development of a Gram-Schmidt projector Pk employing a rank-1 perturbation of the identity matrix. Designed to maximize the high memory bandwidth utilization of the AMD MI-250X GPU, this approach includes a strategy for treating the unit diagonal that minimizes memory reads, leading to a 25\% increase in computational efficiency. The application of perturbation theory further ensures that orthogonality loss is limited to \scrO (\epsilon ) k, where k is the number of iterations. Additionally, we introduce a mixed AMG V-cycle strategy combining ILU(0) and L1-Jacobi smoothers, which achieves a 30--50\% reduction in GPU compute times compared to conventional methods, while maintaining low backward error. This strategy, alongside our novel treatment of the diagonal in triangular matrices, marks a substantial increase in AMG efficicency for GPU systems. We believe that these contributions represent a significant advance in optimizing GMRES+AMG algorithms for GPU computations. The empirical results demonstrate notable speed increments and maintain rigorous backward error bounds, underscoring the potential of our methods to substantially increase computational efficiency in large-scale scientific applications.

One of the limitations of recycled GCRO methods is the large amount of computation required to orthogonalize the basis vectors of the newly generated Krylov subspace for the approximate solution when combined with those of the recycle subspace. Recent advancements in low synchronization Gram-Schmidt and generalized minimal residual algorithms, Swirydowicz et al.~\cite{2020-swirydowicz-nlawa}, Carson et al. \cite{Carson2022}, and Lund \cite{Lund2022}, can be incorporated, thereby mitigating the loss of orthogonality of the basis vectors. An augmented Arnoldi formulation of recycling leads to a matrix decomposition and the associated algorithm can also be viewed as a {\it block} Krylov method. Generalizations of both classical and modified block Gram-Schmidt algorithms have been proposed, Carson et al.~\cite{Carson2022}. Here, an inverse compact $WY$ modified Gram-Schmidt algorithm is applied for the inter-block orthogonalization scheme with a block lower triangular correction matrix $T_k$ at iteration $k$. When combined with a weighted (oblique inner product) projection step, the inverse compact $WY$ scheme leads to significant (over 10$\times$ in certain cases) reductions in the number of solver iterations per linear system. The weight is also interpreted in terms of the angle between restart residuals in LGMRES, as defined by Baker et al.\cite{Baker2005}. In many cases, the recycle subspace eigen-spectrum can substitute for a preconditioner.

Anderson Acceleration (AA) is a method to accelerate the convergence of fixed point iterations for nonlinear, algebraic systems of equations. Due to the requirement of solving a least squares problem at each iteration and a reliance on modified Gram-Schmidt for updating the iteration space, AA requires extra costly synchronization steps for global reductions. Moreover, the number of reductions in each iteration depends on the size of the iteration space. In this work, we introduce three low synchronization orthogonalization algorithms into AA within SUNDIALS that reduce the total number of global reductions per iteration to a constant of 2 or 3, independent of the size of the iteration space. A performance study demonstrates the reduced time required by the new algorithms at large processor counts with CPUs and demonstrates the predicted performance on multi-GPU architectures. Most importantly, we provide convergence and timing data for multiple numerical experiments to demonstrate reliability of the algorithms within AA and improved performance at parallel strong-scaling limits.

The parallel strong-scaling of iterative methods is often determined by the number of global reductions at each iteration. Low-synch Gram–Schmidt algorithms are applied here to the Arnoldi algorithm to reduce the number of global reductions and therefore to improve the parallel strong-scaling of iterative solvers for nonsymmetric matrices such as the GMRES and the Krylov–Schur iterative methods. In the Arnoldi context, the QR factorization is “left-looking” and processes one column at a time. Among the methods for generating an orthogonal basis for the Arnoldi algorithm, the classical Gram–Schmidt algorithm, with reorthogonalization (CGS2) requires three global reductions per iteration. A new variant of CGS2 that requires only one reduction per iteration is presented and applied to the Arnoldi algorithm. Delayed CGS2 (DCGS2) employs the minimum number of global reductions per iteration (one) for a one-column at-a-time algorithm. The main idea behind the new algorithm is to group global reductions by rearranging the order of operations. DCGS2 must be carefully integrated into an Arnoldi expansion or a GMRES solver. Numerical stability experiments assess robustness for Krylov–Schur eigenvalue computations. Performance experiments on the ORNL Summit supercomputer then establish the superiority of DCGS2 over CGS2.

. The GMRES algorithm of Saad and Schultz (1986) for nonsymmetric linear systems relies on the Arnoldi expansion of the Krylov basis. The algorithm computes the QR factorization of the matrix B = [ r 0 ,AV k ] at each iteration. Despite an O ( ε ) κ ( B ) loss of orthogonality, the modified Gram-Schmidt (MGS) formulation was shown to be backward stable in the seminal papers by Paige, et al. (2006) and Paige and Strakoˇs (2002). Classical Gram-Schmidt (CGS) exhibits an O ( ε ) κ 2 ( B ) loss of orthogonality, whereas DCGS-2 (CGS with delayed reorthogonalization) reduces this to O ( ε ) in practice (without a formal proof). We present a post-modern (viz. not classical) GMRES algorithm based on Ruhe (1983) and the low-synch algorithms of ´Swirydowicz et al (2020) that achieves O ( ε ) (cid:107) A v k (cid:107) 2 /h k +1 ,k loss of orthogonality. By projecting the vector A v k , with Gauss-Seidel relaxation, onto the orthogonal complement of the space spanned by the computed Krylov vectors V k where V Tk V k = I + L k + L Tk , we can further demonstrate that the loss of orthogonality is at most O ( ε ) κ ( B ). For a broad class of matrices, unlike MGS-GMRES, significant loss of orthogonality does not occur and the relative residual no longer stagnates for highly non-normal systems. The Krylov vectors remain linearly independent and the smallest singular value of V k is not far from one. We also demonstrate that Henrici’s departure from normality of the lower triangular matrix T k ≈ ( V Tk V k ) − 1 in the modified Gram-Schmidt projector P = I − V k T k V Tk is an appropriate quantity for detecting the loss of orthogonality. Our new algorithm results in an almost symmetric correction matrix T k .

NREL: M.A. Sprague (PI), S. Ananthan, R. Binyahib, M. Brazell, M. Henry de Frahan, R. King, P. Mullowney, J. Rood, A. Sharma, S. Thomas, G. Vijayakumar SNL: P. Crozier, L. Berger-Vergiat, L. Cheung, D. Dement, N. de Velder, D. Glaze, J. Hu, R. Knaus, D. Lee, N. Matula, T. Okusanya, J. Overfelt, S. Rajamanickam, P. Sakievich, T. Smith, J. Vo, A. Williams, I. Yamazaki ORNL: J. Turner, A. Prokopenko, R. Wilson UTA: R. Moser, J. Melvin Parallel Geometric Algorithms: J. Sitaraman

The linear equations that arise in interior methods for constrained optimization are sparse symmetric indefinite and become extremely ill-conditioned as the interior method converges. These linear systems present a challenge for existing solver frameworks based on sparse LU or LDL^T decompositions. We benchmark five well known direct linear solver packages using matrices extracted from power grid optimization problems. The achieved solution accuracy varies greatly among the packages. None of the tested packages delivers significant GPU acceleration for our test cases.

. A low-synchronization MGS-GMRES Krylov solver employing a truncated Neumann series for the inverse compact WY MGS correction matrix T is presented. A corollary to the backward stability result of Paige et al. [1] establishes that T = I − L k is sufficient for convergence of GMRES when (cid:107) L (cid:107) pF = O ( ε p ) κ pF ( B ), p > 1 where B = [ r 0 , AV m ] with condition number κ ( B ). The columns of the strictly lower triangular matrix L are defined by matrix-vector products of Krylov vectors V T 1: k − 2 v k − 1 . The preconditioner is the classical R¨uge-Stuben AMG algorithm with compatible relaxation and inner-outer Gauss-Seidel smoother. This smoother may also be expressed as a truncated Neumann series. Despite the rapid convergence of GMRES-AMG, the cost of an elliptic pressure solver (e.g. for the Navier-Stokes equations), is still substantial. Drop tolerances are applied to the strictly lower triangular matrices arising in the smoother in order to reduce the number of non-zeros and accelerate the time to solution. The number of small matrix elements are found to increase from fine to coarse levels and thus the efficiency gains are greater for large problems with many levels in the V -cycle. The solver is applied to the pressure continuity equation for the incompressible Navier-Stokes equations. The pressure solve time is reduced considerably without a change in the convergence rate.

Block Gram-Schmidt algorithms serve as essential kernels in many scientific computing applications, but for many commonly used variants, a rigorous treatment of their stability properties remains open. This work provides a comprehensive categorization of block Gram-Schmidt algorithms, particularly those used in Krylov subspace methods to build orthonormal bases one block vector at a time. Known stability results are assembled, and new results are summarized or conjectured for important communication-reducing variants. Additionally, new block versions of low-synchronization variants are derived, and their efficacy and stability are demonstrated for a wide range of challenging examples. Numerical examples are computed with a versatile Matlab package hosted at https://github.com/katlund/BlockStab, and scripts for reproducing all results in the paper are provided. Block Gram-Schmidt implementations in popular software packages are discussed, along with a number of open problems. An appendix containing all algorithms type-set in a uniform fashion is provided.