I'm using an iterative subspace algorithm (dsrrit) to obtain the eigenvalues of an eigenvector equation$$-\nabla^2 \mathbf{x} = \lambda\mathbf{x}$$where $\nabla^2$ is the usual Laplacian operator. The algorithm returns a set of eigenvalues $\lambda$ and a subspace $Q$ that spans the eigenvectors of $\lambda$.The literature statesthe programs do not produce a set of eigenvectors corresponding to the eigenvalues computed ... If explicit eigenvectors are desired, they may be obtained by evaluating the eigenvectors of $T$ and applying (2).but I have plotted the vectors in the returned subspace $Q$ and they do indeed match the expected eigenvectors of $\lambda$. Are there particular conditions that guarantee this? (My eigenvectors are real, and my operator is represented by a symmetric matrix)  
Using a subspace iteration method to obtain eigenvalues. Getting eigenvectors too but I don't understand why
eigenvalues
I haven't read it in detail, but from your linked paper (p. 2):SRRIT attempts  to compute a nested sequence of orthonormal bases of $Q_1, Q_2, \ldots, Q_m$. Specifically, if all goes well, the subroutine produces a matrix Q with orthonormal columns having the property that if $|\lambda_i| > |\lambda_{i+1}|$ then $q_1,\ldots,q_i$ span $Q_i$.What this means is that if your matrix has only real eigenvalues, with $|\lambda_1|\geq|\lambda_2|\geq\cdots\geq|\lambda_n|$, then $q_i$ will be an eigenvector for $\lambda_i$, using the fact that the matrix is symmetric and the eigenvectors will be orthogonal to each other.If $\lambda_i, \lambda_{i+1}$ are two complex conjugate eigenvalues, then the corresponding columns of $Q$ wouldn't be the eigenvectors because they would be real.