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.

(possible) solution

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.

  • 1
    $\begingroup$ When you say it's the memory bottleneck, do you mean bandwidth/speed or capacity/size? If it's bandwidth, moving it to disk is going to make things worse. If it's a size thing.. are you sure this is really an issue? Storing a single matrix per element, of fixed size, is not really asymptotically worse than storing the state vector/iterate vector itself (though certainly worse by a large constant factor). If storing them really is using too much memory, I think you would be better served by looking at distributed memory parallelism than out-of-core techniques. $\endgroup$ – rchilton1980 Sep 18 '18 at 14:33
  • $\begingroup$ @rchilton1980 I mean capacity/size. As for the matrix vs. vector, yes, the large constant factor is a problem, and only gets worse with higher polynomial order. I have written a distributed version of the solver that splits the problem between nodes. But on each node, the storage of these inverses or their computation is still the limiting factor. Of course, the distributed approach also introduces a communication overhead between nodes. $\endgroup$ – user3482876 Sep 18 '18 at 14:48
  • 4
    $\begingroup$ It might end up happening that doing the same inverses on the fly is faster than storing them. One keyword you might be looking for regarding the out-of-core solver is io latency hiding. Its efficiency will depend on how your other operations during the iterative solve allow you to incorporate the disk IO without significant slowdown. $\endgroup$ – Anton Menshov Sep 18 '18 at 15:39
  • $\begingroup$ Any repetition in the grid/elements that could be exploited? Can you maybe use a mostly-structured grid? $\endgroup$ – rchilton1980 Sep 18 '18 at 17:04
  • $\begingroup$ @rchilton1980 in general, no; the large problems in which I am interested involve unstructured grids. It's a good suggestion, though. $\endgroup$ – user3482876 Sep 18 '18 at 18:38

There is no faster way to storing things than to put them into memory. If your problem is such that you can't store all of the element matrix inverses in memory, then you have two options: (i) buy more memory for the machines you use, (ii) use more machines. There is also (iii) solve smaller problems. Trying to externalize the matrices to disk will make things vastly slower because disks are at least an order of magnitude slower than memory. Disks are also vastly slower than communication between machines.

But it's worth stepping back. If your solver scheme requires storing lots of local matrix inverses, then maybe that's not a great solver. I'm making the assumption here that in your scheme, storing the inverses is more expensive than storing the local matrices to begin with -- for example, because the local matrices are sparse but the inverse are not. Generally, my experience from 20 years of finite element computations is that if you are running out of memory, then you've already run out of CPU time as well. So I'm a bit surprised that runtimes don't also show up on the list of your worries.

  • 1
    $\begingroup$ just a quick comment on one of the points. If the machine is multicore, I easily run out-of-memory way earlier than I hit the CPU time problem (>1 day) if I don't use distributed cluster but limit myself to a single node. My experience comes from integral equations with dense matrices and the fast methods to accelerate the solutions to underlying problems. $\endgroup$ – Anton Menshov Sep 19 '18 at 23:54
  • $\begingroup$ To your points: the scheme is an HDG scheme, so the element-local contribution matrices to the global linear system require solving a linear system (compare to CG, where elm contributions can simply be computed). Therefore, the element local contributions are also dense. In a matrix-free scheme, these element local contributions must either be stored, or re-computed for each element, for each iteration of an iterative solver. I don't see a way around this, and yes, either time or memory will be limiting depending on the choice above (store vs recompute). $\endgroup$ – user3482876 Sep 20 '18 at 16:47
  • $\begingroup$ Interesting. I've worked with such codes and typically, running out of patience (CPU time) has been the bigger problem. $\endgroup$ – Wolfgang Bangerth Sep 21 '18 at 3:50

One idea would be to store the low-rank approximation of the inverse of the matrix $A^{-1}\approx UV$ using SVD. If $A\in \mathbb{R}^{30\times 30}$, then $U\in \mathbb{R}^{30\times k}\ \& \ V\in \mathcal{R}^{k\times 30} $, where $k\approx 4,5$ is the rank. Use the $UV$ approximation as a preconditioner for solving the system $Ax=b$ using GMRES. You can use vary $k$ as a tuning parameter that trades off memory with the convergence rate.


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.