are needed for the construction of $\boldsymbol{H}$. Therefore, Eq. \[eq:gevp\] is a non-linear problem and must be solved self-consistently. It is possible, and has already been implemented in various codes [@conquest_gillan_2007; @fhiaims_blum_2009; @onetep_skylaris_2005; @siesta_soler_2002], to restrict the computational cost of the construction of $\boldsymbol{H}$ to scale linearly with respect to the system size $N$ for any semi-local and hybrid exchange-correlation functional. In contrast, the solution of a dense eigenproblem (“diagonalization”) scales as $\mathcal{O}(N^3)$, quickly growing to become prohibitive as $N$ increases to large values.

Various strategies exist to facilitate or circumvent the solution of Eq. \[eq:gevp\] [@elsi_yu_2018; @elsi_yu_2019]. (1) The key idea in a conventional dense eigensolver [@numerical_press_2007; @matrix_golub_2013] is tridiagonalization, which brings the original matrix to a tridiagonal form by a series of Householder transformations. This algorithm suffers from the inefficiency of BLAS level-2 matrix-vector operations. New algorithms such as pentadiagonalization [@eigenexa_imamura_2011] and two-stage tridiagonalization [@2stage_lang_1993; @2stage_bischof_1994; @elpa_auckenthaler_2011; @elpa_marek_2014] have been developed, leading to enhanced performance over the conventional one-stage tridiagonalization approach. (2) Iterative eigensolvers [@davidson_davidson_1975; @davidson_sleijpen_1996; @iterative_payne_1992; @iterative_kresse_1996] are commonly employed by DFT codes, particularly those based on plane-wave (PW) basis functions and pseudopotentials. In that case, because of the large number of basis functions, i.e. the large 