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 unknownsThe linear systems which I am solving are symmetric positive definite, but while scipy.sparse doesn'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 numpy array) 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.sparse wraps 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)?
Leveraging scipy for matrix free finite elements
sparse;reference request;scipy;hpc