  We present a non-iterative solver based on the Schur complement method for
sparse linear systems of special form which appear in Quantum Monte-Carlo (QMC)
simulations of strongly interacting fermions on the lattice. While the number
of floating-point operations for this solver scales as the cube of the number
of lattice sites, for practically relevant lattice sizes it is still
significantly faster than iterative solvers such as the Conjugate Gradient
method in the regime of strong inter-fermion interactions, for example, in the
vicinity of quantum phase transitions. The speed-up is even more dramatic for
the solution of multiple linear systems with different right-hand sides. We
present benchmark results for QMC simulations of the tight-binding models on
the hexagonal graphene lattice with on-site (Hubbard) and non-local (Coulomb)
interactions, and demonstrate the potential for further speed-up using GPU.
