dimension of the matrices in Eq. \[eq:gevp\], needed in an accurate calculation, a direct solution of Eq. \[eq:gevp\] is rather infeasible. Iterative eigensolvers are well suited to find a small fraction of low-lying eigenstates, commensurate with the needs of a PW-based code. (3) When using spatially localized basis functions, such as linear combination of atomic orbitals (LCAO), the fraction of needed eigenpairs out of the full matrix dimension can be fairly large. In this scenario, iterative solvers no longer have an advantage over direct eigensolvers. With localized basis functions, locality in the physical system can be translated to sparsity in the $\boldsymbol{H}$ and $\boldsymbol{S}$ matrices. Methods exploiting this sparsity can be formulated as $\mathcal{O}(N) \sim \mathcal{O}(N^2)$ [@linear_goedecker_1999; @linear_bowler_2012; @linear_moussa_2019; @feast_polizzi_2009; @pexsi_lin_2013; @ntpoly_dawson_2018] by circumventing the explicit solution of Eq. \[eq:gevp\]. In particular, linear scaling algorithms in a density matrix formalism have been successfully applied to simulations of one million atoms [@million_bowler_2010; @million_vandevondele_2012]. Despite the success in extreme-scale simulations, reduced scaling methods come with a computational prefactor that is much larger than that of the $\mathcal{O}(N^3)$ diagonalization method. Moreover, the applicability and optimal performance of reduced scaling methods are often limited to some certain problem types. As of today, dense eigensolvers, 