  We propose algorithms for efficient time integration of large systems of
oscillatory second order ordinary differential equations (ODEs) whose solution
can be expressed in terms of trigonometric matrix functions. Our algorithms are
based on a residual notion for second order ODEs, which allows to extend the
``residual-time restarting'' Krylov subspace framework -- which was recently
introduced for exponential and $\varphi$-functions occurring in time
integration of first order ODEs -- to our setting. We then show that the
computational cost can be further reduced in many cases by using our restarting
in the Gautschi cosine scheme. We analyze residual convergence in terms of
Faber and Chebyshev series and supplement these theoretical results by
numerical experiments illustrating the efficiency of the proposed methods.
