An algorithm to compute the extremal eigenvalues of a non-Hermitian matrix, based on Arnoldi-orthogonalization, is employed in the linear stability analysis of the viscoelastic Poiseuille flow at high Reynolds numbers. It is shown that this algorithm is both computationally efficient and accurate in reproducing the most unstable modes in the pseudo-spectrally discretized eigenspectrum of the original problem. The Upper Convected Maxwell (UCM), Oldroyd-B and Chilcott-Rallison models are considered for the linear stability analysis of the high Reynolds number viscoelastic Poiseuille flow. Results for the UCM model show a large destabilization of the flow compared to the Newtonian limit, even for low values of flow elasticity, epsilon = We/Re, of the order 10(-3), realized at high Reynolds numbers for 0(1) We. These results agree qualitatively with those reported by Porteus and Denn, although some quantitative differences exist, especially at high values of elasticity. Furthermore, it is shown that the number of spectral modes necessary to obtain converged results increases substantially as the flow elasticity is increased. A comparison of the linear stability characteristics of the Oldroyd-B and the UCM models has revealed that the presence of a non-zero solvent viscosity has a pronounced stabilizing effect on the flow. Further stabilization occurs through the introduction of a finite molecular extensibility in the Chilcott-Rallison model.