  Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) is widely
used to compute eigenvalues of large sparse symmetric matrices. The algorithm
can suffer from numerical instability if it is not implemented with care. This
is especially problematic when the number of eigenpairs to be computed is
relatively large. In this paper we propose an improved basis selection strategy
based on earlier work by Hetmaniuk and Lehoucq as well as a robust convergence
criterion which is backward stable to enhance the robustness. We also suggest
several algorithmic optimizations that improve performance of practical LOBPCG
implementations. Numerical examples confirm that our approach consistently and
significantly outperforms previous competing approaches in both stability and
speed.
