Appendix A

The pseudo-code for solving the generalized eigenvalue problem $H x = \varepsilon S x$ based on the block update procedure is as follows, where $H$, $S$ and $T$ denote the $N \times N$ Hamiltonian, overlap and kinetic energy matrices respectively. Let $X^{(k)} = \left[ X^{(k)}_1,X^{(k)}_2,\ldots,X^{(k)}_M \right]$ be an $N \times M$ matrix where $X^{(k)}_i$ is the $i$-th column of $X^{(k)}$, and $k$ labels the iteration.

The procedure is then:

$k := 1$ 


Choose $X^{(1)}$ 


Choose convergence tolerance $\epsilon$ 


Orthonormalize $X^{(1)}$ (Gram-Schmidt) 


Calculate $\Omega^{(1)}$ 


Choose $\Omega^{(0)}$ such that $\left\vert \Omega^{(1)} - \Omega^{(0)}\right\vert > \epsilon$ 

$F^{(0)} := 0$ 

$A^{(0)} := 0$ 

$\gamma^{(0)}_1 := 1$ 

 while $\left\vert \Omega^{(k)} - \Omega^{(k-1)} \right\vert > \epsilon$ do 

		 $Y^{(k)} := H X^{(k)}$ 

		 $Z^{(k)} := S X^{(k)}$ 

		 $P^{(k)} := (X^{(k)})^T Y^{(k)}$ 

		 $F^{(k)} := Y^{(k)} - Z^{(k)}P^{(k)}$ 


Solve $(S+T/\tau)B^{(k)} = F^{(k)}$ 

		 $C^{(k)} := (Z^{(k)})^T B^{(k)}$ 

		 $G^{(k)} := B^{(k)} - X^{(k)}C^{(k)}$ 

		 $\gamma^{(k)}_1 := {\mathrm{tr}}\left( (G^{(k)})^T F^{(k)} \right)$

		 $\gamma^{(k)}_2 := {\mathrm{tr}}\left( (G^{(k)})^T F^{(k-1)}\right)$ 

		 $\gamma^{(k)} := \left( \gamma^{(k)}_1 - \gamma^{(k)}_2 \right) /\gamma^{(k-1)}_1$ 

		 $A^{(k)} := -G^{(k)} + \gamma^{(k)} A^{(k-1)}$ 

		 $E^{(k)} := (Z^{(k)})^T A^{(k)}$ 

		 $D^{(k)} := A^{(k)} - X^{(k)}E^{(k)}$ 

		 $\lambda_{\mathrm{opt}} := {\mathrm{linmin}}(X^{(k)},D^{(k)})$ (see Appendix B) 

		 $X^{(k+1)} := X^{(k)} + \lambda_{\mathrm{opt}} D^{(k)}$ 


Orthonormalize $X^{(k+1)}$ (Gram-Schmidt) 


Calculate $\Omega^{(k+1)}$ 

		 $k := k + 1$ 

 end

In some applications, individual eigenvalues are needed. They can be obtained by a subspace rotation method where we simply need to diagonalize the $M \times M$ matrix $X^T H X$. If $U$ diagonalizes $X^T H X$ such that $U^T ( X^T H X) U =$ diag( $\varepsilon_1,\varepsilon_2,\ldots,\varepsilon_M$), we obtain the individual eigenvalues $\{ \varepsilon_i; i = 1,\ldots, M \}$ with corresponding eigenvectors $X' = XU$.