I am trying to solve $A X = B$ for $X$ using some iterative method, where $A$ is a large, square, full rank n-by-n matrix and $X, B$ are also both n-by-n matrices. I'm aware of the large literature on numerical solutions to the case where $X$ and $B$ are both vectors, but haven't been able to find guidance on the more general case.
I could transform the general equation into a linear system, either by solving for each column of $X$ separately or equivalently by solving the vectorized version: $$(I_n \otimes A) x^*= b^*,$$ where $x^*$ and $b^*$ are vectors constructed by stacking the columns of $X$ and $B$. Unfortunately I have about a million columns and each column takes five minutes to solve, and thus solving the system in such a way would take roughly 9 years and 7 months. Is a quicker method available?
You should investigate these keywords: "block Krylov methods" and "Krylov subspace recycling".
Block Krylov methods are designed for your situation, i.e., multiple right-hand sides.
In general, you can use the Krylov subspace generated during the solution of one system to accele-rate the solution of another system with the same coefficient matrix. This procedure is known as Krylov subspace recycling.