I am interested in solving many linear systems $Ax = b$, where $A$ is symmetric positive definite and small (i.e. less than 25,000 rows) --- $b$ will be changing. We can assume that $A$ arises from the FEM discretization of the elasticity equation in 3D on an unstructured mesh for an arbitrarily-shaped domain.
I am looking for a method that uses less memory than a sparse solver and such that one solution takes less than twice the time for a forward-backward.
Which algorithm would you recommend under these assumptions?
I would be grateful if you could mention why it would be efficient for small matrices (i.e. less than 25,000 rows). For such matrix sizes, I believe that, even if the algorithm is $O(n)$, the constant in front of $n$ matters a lot.