Given a real symmetric matrix \(\mathbf{A} = {\{a_{ij}\}}_{n \times n}\) of order \(n\), calculates the \(k\) largest eigenvalues, \(\lambda _1 \ge \lambda _2 \ge ,..., \ge \lambda _k\) and their associated eigenvectors \(\mathbf{x}_i , i=1,2,..., k\), where \(\mathbf{Ax} _i = \lambda _i \mathbf{x} _i\). The subroutine uses the method of simultaneous iteration and also allows the user to take advantage of sparsity. Eigensolutions from other parts of the spectrum can be obtained by using shifts of the form \(\mathbf{A} - \beta \mathbf{I}\) optionally combined with inversion.