known as strong interaction, which describes the interactions among quarks, one of the constituents of matter. Lattice QCD is the tool for non-perturbative numerical calculations of these interactions on a Euclidean space-time lattice [@Wilson:1974sk]. The heart of the computations is the solution of the lattice-Dirac equation, which translates to a linear system of equations $Mx = b$, often for a large number of right hand sides [@Foley:2005ac]. The Dirac operator $M$ is $\gamma_5$-Hermitian, or $\gamma_5 M = M^H \gamma_5$, where, in one representation, $\gamma_5$ is a diagonal matrix with 1 and -1 on the diagonal. Also, $M = m_q I- D$, where $m_{q}$ is a parameter related to the quark mass and $D$ is an operator. In addition to solving linear systems of equations, many current approaches [@Foley:2005ac; @Edwards:2001ei; @DeGrand:2004qw; @Giusti:2004yp] require the solution of the eigenvalue problem $(\gamma_5 M)u_i = \lambda_i u_i$, for 100-200 smallest magnitude eigenvalues, $\lambda_i$, and their eigenvectors, $u_i$ (together we call them eigenpairs). Beyond the very large dimension and number of right hand sides, $M$ becomes increasingly ill-conditioned as $m_{q} \rightarrow m_{critical}$. In lattice QCD this is known as critical slowdown and is a limiting computational factor.

Traditionally these linear systems are solved by applying 