  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.
