Resolvent analysis has found applications in several areas of fluid mechanics, providing physical insight into both laminar and turbulent flows. In spite of such fact, the global (3D) resolvent is computationally expensive, which limits the size of the domain and the Reynolds number of the flows which can be considered. In this work, we derive a parabolic resolvent approach, which enables a significant increase in the computational efficiency of the calculation, for streaky structures in boundary layer flows. The computational speedup depends on the size of the problem and could be of more than one order of magnitude for the same accuracy as the global calculation. The method is derived based on an optimiza-tion method via the Lagrange multipliers over the linearized boundary layer equations and it is coupled to a Krylov-Arnoldi decomposition to the computation of suboptimals. The application of the method is exemplified for two problems: a Falkner-Skan boundary layer, where we obtain trends for both the optimals and suboptimals, and a turbulent boundary layer, where characteristics such as the double peak in the spectrum and the characteristic inner and outer length scales can be recovered when a variable eddy viscosity is considered. In both cases, a scaling is found for the dominant gain, given in terms of the fourth power of the Reynolds number, defined in terms of the relevant scale for the problem, the displacement thickness, and the modified Rotta-Clauser parameter for the laminar and turbulent boundary layers, respectively. For the laminar case, we further demonstrate that a forcing limited to the free-stream region is capable of generating streaky structures inside the boundary layer, a relevant feature for free-stream turbulence-induced transition.