  Recently a new algorithm for model reduction of second order linear dynamical
systems with proportional damping, the Adaptive Iterative Rational Global
Arnoldi (AIRGA) algorithm (Bonin et. al., 2016), has been proposed. The main
computational cost of the AIRGA algorithm is solving a sequence of linear
systems. Usually, direct methods (e.g., LU) are used for solving these systems.
As model sizes grow, direct methods become prohibitively expensive. Iterative
methods (e.g., Krylov) scale well with size, and hence, are a good choice with
an appropriate preconditioner.
  Preconditioned iterative methods introduce errors in linear solves because
they are not exact. They solve linear systems up to a certain tolerance. We
prove that, under mild conditions, the AIRGA algorithm is backward stable with
respect to the errors introduced by these inexact linear solves. Our first
assumption is use of a Ritz-Galerkin based solver that satisfies few extra
orthogonality conditions. Since Conjugate Gradient (CG) is the most popular
method based upon the Ritz-Galerkin theory, we use it. We show how to modify CG
to achieve these extra orthogonalities.
  Modifying CG with the suggested changes is non-trivial. Hence, we demonstrate
that using Recycling CG (RCG) helps us achieve these orthogonalities with no
code changes. The extra cost of orthogonalizations is often offset by savings
because of recycling. Our second and third assumptions involve existence,
invertibility and boundedness of two matrices, which are easy to satisfy.
  While satisfying the backward stability assumptions, by numerical experiments
we show that as we iteratively solve the linear systems arising in the AIRGA
algorithm more accurately, we obtain a more accurate reduced system. We also
show that recycling Krylov subspaces helps satisfy the backward stability
assumptions (extra-orthogonalities) at almost no extra cost.
