Efficient iterative refinement and regularized LDL solvers with Multiple Right-Hand Sides

[poulson]

Current iterative refinement and regularized solvers, e.g., reg_ldl::SolveAfter, individually solve for each of the right-hand sides or simply only support a single right-hand side. Given that Homogeneous Self-Dual embeddings of linear conic programs introduce an extra linear solve, which can happen simultaneously with another solve, it makes sense to avoid loading the sparse matrix data through the memory hierarchy twice as many times as necessary.

It might be worth having specialized implementations for one, two, and n right-hand sides to efficiently handle deflation.