& = \boldsymbol{P} \boldsymbol{\tilde{H}} \boldsymbol{P}^* , \label{eq:2stage_fwd1} \\
\boldsymbol{T} & = \boldsymbol{Q} \boldsymbol{B} \boldsymbol{Q}^* , \label{eq:2stage_fwd2} \\
\boldsymbol{T} \boldsymbol{X} & = \boldsymbol{X} \boldsymbol{\Sigma} , \label{eq:2stage_solve} \\
\boldsymbol{Y} & = \boldsymbol{Q}^* \boldsymbol{X} , \label{eq:2stage_bkwd1} \\
\boldsymbol{\tilde{C}} & = \boldsymbol{P}^* \boldsymbol{Y} . \label{eq:2stage_bkwd2}\end{aligned}$$

![image](fig/2stage.png){width="99.00000%"}

This two-stage tridiagonalization approach is implemented in several linear algebra software packages [@elpa_auckenthaler_2011; @elpa_marek_2014; @magma_tomov_2010; @lapack_anderson_1999; @mkl_arturov_2018], including a high-performance distributed-memory implementation in the ELPA library [@elpa_auckenthaler_2011; @elpa_marek_2014]. The introduction of the banded matrix stage leads to faster computation compared to the one-stage method, for the transformation in Eq. \[eq:2stage\_fwd1\] mostly involves highly efficient BLAS level-3 matrix-matrix operations, and the transformation in Eq. \[eq:2stage\_fwd2\] only works on a sparse banded matrix $\boldsymbol{B}$ instead of a full matrix. The solution of Eq. \[eq:2stage\_solve\] is accelerated in ELPA by extending the divide-and-conquer symmetric tridiagonal eigensolver [@dc_cuppen_1980; @dc_gu_1995; @dc_tisseur_1999] such that unwanted eigenvectors are not computed [@elpa_auckenthaler_2011; @elpa_auckenthaler_2013]. Regarding the back-transformation of eigenvectors, it is rather difficult to directly use BLAS routines for Eq. \[eq:2stage\_bkwd1\]. Manually optimized “kernels” written in architecture-specific instruction sets are employed for this particular step [@elpa_auckenthaler_2011; @elpa_kus_2019b; @elpa_auckenthaler_2013]. The ELPA2 solver is highly scalable on massively parallel, distributed-memory architectures. It avoids global MPI communications by using a 