Rational hybrid Monte Carlo with block solvers and multiple pseudofermions

The dominant cost of most lattice QCD simulations is the inversion of the Dirac operator required to calculate the force term in the rational hybrid Monte Carlo (RHMC) update. One way to improve this situation is to use multiple pseudofermions, which reduces the size and variance of this force and hence allows a larger integration step size to be used. This means fewer force term calculations are required, but at the cost of having to invert the Dirac operator for each pseudofermion field. This bottleneck can be addressed: recently there has been renewed interest in the use of block Krylov solvers, which can solve multiple right-hand-side vectors with significantly fewer iterations than are required if each vector is solved using a separate Krylov solver. We combine these two ideas, achieving a significant speed-up of RHMC lattice QCD simulations.