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

  • $\begingroup$ I know trilinos has python wrappers that work well to supplement scipy/numpy $\endgroup$
    – KyleW
    Jul 25, 2017 at 15:44

2 Answers 2


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.


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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