Generalized perturbation theory (GPT) has been recognized as the most computationally efficient approach for performing sensitivity analysis for models with many input parameters, which renders forward sensitivity analysis computationally overwhelming. In critical systems, GPT involves the solution of the adjoint form of the eigenvalue problem with a response-dependent fixed source. Although conceptually simple to implement, most neutronics codes that can solve the adjoint eigenvalue problem do not have a GPT capability unless envisioned during code development. We introduce in this manuscript a reduced-order modeling approach based on subspace methods that requires the solution of the fundamental adjoint equations but allows the generation of response sensitivities without the need to set up GPT equations, and that provides an estimate of the error resulting from the reduction. Moreover, the new approach solves the eigenvalue problem independently of the number or type of responses. This allows for an efficient computation of sensitivities when many responses are required. This paper introduces the theory and implementation details of the GPT-free approach and describes how the errors could be estimated as part of the analysis. The applicability is demonstrated by estimating the variations in the flux distribution everywhere in the phase space of a fast critical sphere and a high-temperature gas-cooled reactor prismatic lattice. The variations generated by the GPT-free approach are benchmarked to the exact variations generated by direct forward perturbations.