In the case of a finite element code, I have many small (order of 30x30) matrix inverses (or LU factorizations), one per finite element. These matrix inverses never change and must be applied to vectors repeatedly in the context of a matrix-free solver.
Storing all these inverses is the memory bottleneck of the code, because they have to be applied for each finite element in the mesh. On the other hand, computing the inverses 'on the fly' (using LAPACK or Numpy) is not practical, because these inverses must be applied at each iteration of a sparse iterative solve, so each time-step would require
(number of CG iterations) * (number of elements) * (time of inverse)
which is too expensive, even for problems far smaller than those I am interested in.
I am considering performing the matrix inversions or LU factorizations once, and then storing them to disk, rather than keeping them in memory. Then the idea is to read the inverses into memory either one at a time or by chunks, and apply them, avoiding the repeated inversions. But I know that reading from disk can be expensive.
Are there good or accepted ways of doing this? Most out-of-core papers I've looked at are considering direct solves or much more complicated problems, so I would appreciate any pointers to relevant references as well.
EDIT: The code is a mix of Python and C, and I've found some posts indicating that perhaps the HDF5 library is a good alternative to rolling my own cache, as suggested above.