The traditional approach to investigating the stability of a physical system is to linearize the equations about a steady base solution, and to examine the eigenvalues of the linearized operator. Over the past several decades, it has been recognized that this approach only determines the asymptotic stability of the system, and neglects the possibility of transient perturbation growth arising due to the nonnormality of the system. This observation motivated the development of a more powerful generalized stability theory (GST), which focuses instead on the singular value decomposition (SVD) of the linearized propagator of the system. While GST has had significant successes in understanding the stability of phenomena in geophysical fluid dynamics, its more widespread applicability has been hampered by the fact that computing the SVD requires both the tangent linear operator and its adjoint: deriving the tangent linear and adjoint models is usually a considerable challenge, and manually embedding them inside an eigensolver is laborious. In this paper, we present a framework for the automation of generalized stability theory, which overcomes these difficulties. Given a compact high-level symbolic representation of a finite element discretization implemented in the FEniCS system, efficient C++ code is automatically generated to assemble the forward, tangent linear, and adjoint models; these models are then used to calculate the optimally growing perturbations to the forward model, as well as their growth rates. By automating the stability computations, we hope to make these powerful tools a more routine part of computational analysis. The efficiency and generality of the framework are demonstrated, with applications drawn from geophysical fluid dynamics, phase separation, and quantum mechanics.