[SciPy-user] scipy.sparse.linalg.cg not thread safe?

[Nathan Bell]

On Tue, Jan 27, 2009 at 4:31 PM, Sturla Molden <sturla at molden.no> wrote:
>
> How would this be performance wise? Some iterative methods are fast in
> Python, others are not. Why not protect Fortran code unsafe for threads
> with a global lock in the Python module? It would not be any worse than
> the GIL, which would affect pure Python code.
>

I don't have an opinion on the locking issue, but the dominant cost in
most iterative methods for linear systems is the cost of the sparse
matrix-vector products (for y = A*x for sparse A).  A smaller amount
of time is spent in level 1 BLAS operations like axpy() and norm().
All of these map efficiently to existing Python + SciPy functionality,
so there's little overhead.

IMO the advantage of the pure Python approach is also evident.
Compare the Python code that interfaces to the Fortran CG
implementation to the pure Python + SciPy implementation of the
*entire* algorithm:
http://projects.scipy.org/scipy/scipy/browser/trunk/scipy/sparse/linalg/isolve/iterative.py#L196
http://code.google.com/p/pyamg/source/browse/trunk/pyamg/krylov/cg.py#71

I'm definitely in favor of parallelizing as much as of SciPy as
possible.  Now that most compilers support OpenMP it would be fairly
straightforward to parallelize the C++ code that implements sparse
matrix-vector multiplication (among other things).  In conjunction, we
ought to add the necessary compiler flags to setuptools so that
OpenMP-enabled sources are handled correctly.

-- 
Nathan Bell wnbell at gmail.com
http://graphics.cs.uiuc.edu/~wnbell/