  One of the most computationally expensive steps of the low-rank ADI method
for large-scale Lyapunov equations is the solution of a shifted linear system
at each iteration. We propose the use of the extended Krylov subspace method
for this task. In particular, we illustrate how a single approximation space
can be constructed to solve all the shifted linear systems needed to achieve a
prescribed accuracy in terms of Lyapunov residual norm. Moreover, we show how
to fully merge the two iterative procedures in order to obtain a novel,
efficient implementation of the low-rank ADI method, for an important class of
equations. Many state-of-the-art algorithms for the shift computation can be
easily incorporated into our new scheme, as well. Several numerical results
illustrate the potential of our novel procedure when compared to an
implementation of the low-rank ADI method based on sparse direct solvers for
the shifted linear systems.
