  We present the Alternating Anderson-Richardson (AAR) method: an efficient and
scalable alternative to preconditioned Krylov solvers for the solution of
large, sparse linear systems on high performance computing platforms.
Specifically, we generalize the recently proposed Alternating Anderson-Jacobi
(AAJ) method (Pratapa et al., J. Comput. Phys. (2016), 306, 43--54) to include
preconditioning, discuss efficient parallel implementation, and provide serial
MATLAB and parallel C/C++ implementations. In serial applications to
nonsymmetric systems, we find that AAR is comparably robust to GMRES, using the
same preconditioning, while often outperforming it in time to solution; and
find AAR to be more robust than Bi-CGSTAB for the problems considered. In
parallel applications to the Helmholtz and Poisson equations, we find that AAR
shows superior strong and weak scaling to GMRES, Bi-CGSTAB, and Conjugate
Gradient (CG) methods, using the same preconditioning, with consistently
shorter times to solution at larger processor counts. Finally, in massively
parallel applications to the Poisson equation, on up to 110,592 processors, we
find that AAR shows superior strong and weak scaling to CG, with shorter
minimum time to solution. We thus find that AAR offers a robust and efficient
alternative to current state-of-the-art solvers, with increasing advantages as
the number of processors grows.
