Linear stability analysis of incompressible viscoelastic flows based on normal mode expansions of the eigenfunctions requires the numerical solution of a generalized eigenvalue problem (GEVP). The complex boundary layer structure of the leading eigenfunctions and the singular character of the continuous set of eigenvalues, necessitate the use of fine mesh sizes, leading to large algebraic GEVPs. In this paper, we present a submatrix-based transformation of the linearized equations (SubTLE) that converts the GEVP into a simple eigenvalue problem (EVP) of half the original dimension for the purely elastic isothermal and non-isothermal flows of an Oldroyd-B liquid. This leads to significant (up to an order of magnitude) reduction in the CPU time and memory required for the solution of the EVP. This is illustrated in the context of isothermal and non-isothermal shear flows.