The parabolized stability equation (PSE) is a widely used, efficient method to calculate the evolution of streamwise traveling instabilities in spatially developing, weakly nonparallel flows. Although the PSE method is very economical, computational time can be significant due to the repeated solution of linear systems of equations in each marching step. The linear systems of equations are typically solved by LU factorization due to the moderate size and the sparsity of the coefficient matrices. In this paper, instead of repeated calculation of the LU factorization, a single LU factorization is calculated in each marching step, and used as a preconditioner for an iterative solver. Numerical experiments are conducted on the incompressible PSE, nonlinear PSE (NLPSE) with explicit discretization of the nonlinear terms, and NLPSE with implicit treatment of the nonlinearities. It is shown that the time spent on the solution of the linear systems of equations was reduced by a factor of 2.5, 4.5-7.5, and 20-60 for the three cases, respectively. The methodology can be generalized to the compressible PSE equations, and extended to similar boundary layer stability problems, or slowly varying parabolic partial differential equations.