I have a large number of systems of the form:
To solve for a large numbers of such $b_i\;1\leq i \leq k$ but where $A$ is fixed (A is a rank $p$ general --i.e. non sparse, non PSD-- matrix).
I can solve them individually using an LU decomposition (costs $O(p^3)$) but was wondering whether there is a more efficient way to get all of the $k$ vector of solution $x^*_i$ than solving these $k$ systems independently?
The typical range of values of $p$, $k$ i'm considering are in the 10-100 (typically, i'm expecting $k\approx p$).
A pointer to a c++ implementation of whatever method is proposed would also be greatly appreciated.