This paper describes the development and initial application of an adjoint harmonic balance (HB) solver. The HB method is a numerical method formulated in the frequency domain which is particularly suitable for the simulation of periodic unsteady flow phenomena in turbomachinery. Successful applications of this method include unsteady aerodynamics as well as aeroacoustics and aeroelasticity. Here, we focus on forced response due to the interaction of neighboring blade rows. In the simulation-based design and optimization of turbomachinery components, it is often helpful to be able to compute not only the objective values—e.g., performance data of a component—themselves but also their sensitivities with respect to variations of the geometry. An efficient way to compute such sensitivities for a large number of geometric changes is the application of the adjoint method. While this is frequently used in the context of steady computational fluid dynamics (CFD), it becomes prohibitively expensive for unsteady simulations in the time domain. For unsteady methods in the frequency domain, the use of adjoint solvers is feasible but still challenging. The present approach employs the reverse mode of algorithmic differentiation (AD) to construct a discrete adjoint of an existing HB solver in the framework of an industrially applied CFD code. The paper discusses implementational issues as well as the performance of the adjoint solver, in particular regarding memory requirements. The presented method is applied to compute the sensitivities of aeroelastic objectives with respect to geometric changes in a turbine stage.