(compared to the original problem) solved. Then, eigenvectors of $\boldsymbol{T}$ are back-transformed to obtain eigenvectors of $\boldsymbol{\tilde{H}}$. This algorithm is adopted by a variety of dense linear algebra packages, such as LAPACK [@lapack_anderson_1999] targeting sequential and shared-memory parallel architectures; cuSOLVER [@cusolver] and MAGMA [@magma_tomov_2010; @magma_dongarra_2014] targeting shared-memory architectures with GPU accelerators; ScaLAPACK [@scalapack_blackford_1997] and Elemental [@elemental_poulson_2013] targeting distributed-memory parallel architectures. However, as mentioned above, the tridiagonalization step makes extensive use of memory-bound, BLAS level-2 matrix-vector operations, whose performance is limited on modern computer architectures.

The two-stage tridiagonalization algorithm proposed by Bischof, Sun, and Lang [@2stage_lang_1993; @2stage_bischof_1994] is an established alternative to the conventional one-stage method. As shown in Eq. \[eq:2stage\] below and further illustrated in Fig. \[fig:2stage\], the tridiagonalization of the full matrix $\boldsymbol{\tilde{H}}$ is carried out in two transformations. The first transformation $\boldsymbol{P}$ reduces $\boldsymbol{\tilde{H}}$ to a banded matrix $\boldsymbol{B}$, and the second transformation $\boldsymbol{Q}$ reduces $\boldsymbol{B}$ to a tridiagonal matrix $\boldsymbol{T}$. The eigenvalues $\boldsymbol{\Sigma}$ and eigenvectors $\boldsymbol{X}$ of $\boldsymbol{T}$ are solved as done in the one-stage method. The back-transformation of eigenvectors is also carried out in two steps. $\boldsymbol{X}$ is first back-transformed to $\boldsymbol{Y}$, the eigenvectors of $\boldsymbol{B}$, then to $\boldsymbol{\tilde{C}}$, the eigenvectors of $\boldsymbol{\tilde{H}}$.

\[eq:2stage\] $$\begin{aligned}
\boldsymbol{B} 