Leveraging scipy for matrix free finite elements

Question

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.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)?

Answer 1

I would say that the implementation + verification + unit testing would take you more than just 3 weeks. Although, if you are planning to invest that time, you might add that capabilities to scipy.sparse or scikits-sparse.

Regarding symmetric sparse matrices, you can check Pysparse. It has Sparse Skyline format (SSS) that is used for symmetric matrices. You can also directly use the Low-Level Sparse Matrices that they provide.

This package is used by SfePy that is a FEM package for Python.

Answer 2

It turns out that scipy does indeed support this type of overloading.

One simply needs to write a class inheriting from scipy.sparse.linalg.LinearOperator and implement the matvec method.

Doing so allows use of all the scipy.sparse linear solvers, but with a matrix vector routine that can be implemented as the user desires (in my case, which applies each contribution independently and takes advantage of symmetry).