For my tasks, I had some success with algebraic multigrid solvers as preconditioner, for example from AMGCL or PyAMG. They are also reasonably easy to get start

For my tasks, I had some success with algebraic multigrid solvers as preconditioner, for example from AMGCL or PyAMG. They are also reasonably easy to get started with. https://github.com/pyamg/pyamg https://github.com/pyamg/pyamg https://github.com/ddemidov/amgcl https://github.com/ddemidov/amgcl But I only have to deal with positive definite systems, so YMMV. I am not sure whether those libraries can deal with multiple right-hand sides, but most complexity is in the preconditioners anyway.