  This study computes the gradient of a function of numerical solutions of
ordinary differential equations (ODEs) with respect to the initial condition.
The adjoint method computes the gradient approximately by solving the
corresponding adjoint system numerically. In this context, Sanz-Serna [SIAM
Rev., 58 (2016), pp. 3--33] showed that when the initial value problem is
solved by a Runge--Kutta (RK) method, the gradient can be exactly computed by
applying an appropriate RK method to the adjoint system. Focusing on the case
where the initial value problem is solved by a partitioned RK (PRK) method,
this paper presents a numerical method, which can be seen as a generalization
of PRK methods, for the adjoint system that gives the exact gradient.
