Solving the linear quadratic regulator (LQR) problem for partial differential equations (PDEs) leads to many computational challenges. The primary challenge comes from the fact that discretization methods for PDEs typically lead to very large systems of differential or differential algebraic equations. These systems are used to form algebraic Riccati equations involving high rank matrices. Although we restrict our attention to control problems with small numbers of control inputs, we allow for potentially high order control outputs. Problems with this structure appear in a number of practical applications yet no suitable algorithm exists. We propose and analyze solution strategies based on applying model order reduction methods to Chandrasekhar equations, Lyapunov/Sylvester equations, or combinations of these equations. Our numerical examples illustrate improvements in computational time up to several orders of magnitude over standard tools (when these tools can be used). We also present examples that cannot be solved using existing methods. These cases are motivated by flow control problems that are solved by computing feedback controllers for the linearized system.