I will present a method for calculating eigenvectors of the staggered Dirac operator based on the Golub-Kahan-Lanczos bidiagonalization algorithm. Instead of using orthogonalization during the bidiagonalization procedure to increase stability, we choose to stabilize the method by combining it with an outer iteration that refines the approximate eigenvectors obtained from the inner bidiagonalization procedure. I will discuss the performance of the current implementation using QEX and compare with other methods.