Say, we have a Hamiltonian which for simplicity does not mix particle hole sectors. It is just a simple Hamiltonian in real space as shown,
$H=\sum_{ij,\sigma} A(i,j)(c_{i\sigma}^{\dagger}c_{j\sigma} +h.c.) $
In a specific spin sector (as spin sectors are decoupled) , we use the transformation $c_{i\sigma}=\sum_{n} \alpha_{n\sigma}(i)d_{n\sigma}$ to diagonalize the Hamiltonian to $H_{d}=\sum_{n,\sigma}E_{n\sigma}d_{n\sigma}^{\dagger}d_{n\sigma}$.
If we use LAPACK subroutine (dsyev) to diagonalize the Hamiltonian we will get a matrix consisting of the eigenvectors and corresponding eigenvalues. LAPACK organizes eigenvalues and corresponding vectors in increasing order of their value. Now is there any way to label the eigenvalues. Let me be specific. So in a mean field calculation where you solve equations self consistently, you have to refer to the eigenvalues $E_{n\sigma}$ in the self consistent equations (which you find out analytically). Now in each iteration LAPACK mixes the order of eigenvalues, so we can not refer to the first eigenvalue that it returns in each iteration as $E_{1}$.
Is there a way out?