This will be a very general question.
I have a 3D finite element code in Python which I would like to extend to handle "large" problems (~10^8 unknowns in the global system). Right now I am using the
scipy.sparse library, which gives decent performance for iterative solvers, but I'm finding the following problems:
- I'm quickly running into memory limitations for problems larger than 10^6 unknowns
- The linear systems which I am solving are symmetric positive definite, but while
scipy.sparsedoesn't seem to have a storage format which knows about symmetry, so I think I am storing many more entries than necessary.
- Even before solving the global system, computing the element-local contributions (held in a multidimensional
numpyarray) are abutting the memory requirements.
Therefore, it seems clear to me that I need to either write my own matrix-free iterative solver, or to use an external library with that functionality. My questions:
- I understand that
scipy.sparsewraps lower level routines such as lapack. What I need to write basically overloads sparse matrix vector multiplication. If I write this in Python (or maybe in C and wrap it with Python), do I have any hope of attaining decent performance, or is it simply necessary to be able to call these lower level libraries? I have no intuition for this, and don't want to spend 3 weeks writing my own matrix-free A*x routine which is too slow.
- Is it a good idea to write this routine in Python? Or do I need to write it in C and wrap it?
- Do there exist libraries which support this functionality? It seems unlikely, since a code-specific knowledge of how to compute A*x without assembling A would be necessary.
- Rather than write a matrix-free routine, could I write an out-of-core solver? Do such approaches scale to problems like mine?
- I have access to a cluster with many multi-core compute nodes. Eventually I would like to parallelize this implementation. Are there good tools or references for implementation guidelines?
Or is there a good way to handle the above while using only
scipy? For example, is it possible to provide
scipy.sparse.linalg.cg a pointer to a function which computes A * x (obviously the cg routine is limited by the A*x speed, but it would be nice to avoid coding CG, GMRES, etc. myself)?