[SciPy-Dev] Gmres with Sparse Complex Matrices gives Segmentation fault (or NULL result)

[Martin Campos Pinto]

Hi all, I am having some troubles using iterative solvers with sparse complex matrices. (I first posted that issue on the SciPy-User mailing list but I guess this one is more appropriate.)

So, here is an example with gmres: the script

from scipy import sparse
from numpy.linalg import norm
from numpy.random import rand
import scipy.sparse.linalg as spla

C = sparse.lil_matrix((10, 10), dtype=complex)
C.setdiag(rand(10))
C[0,3] = 1j
C = C.tocsr()
c = rand(10)+rand(10)*1j

x = spla.spsolve(C, c)
print "spsolve: norm(C*x - c) = ", norm(C*x - c)

(x,info) = spla.gmres(C, c)
print "gmres: norm(C*x - c) = ", norm(C*x - c)

gives as output:

spsolve: norm(C*x - c) =  1.57621370006e-16
Segmentation fault


Actually, sometimes I get this error message instead of a Segmentation fault:

Traceback (most recent call last):
  File "test_gmres_cplx.py", line 45, in <module>
    (x,info) = spla.gmres(C, c)
  File "/Library/Frameworks/EPD64.framework/Versions/7.2/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/iterative.py", line 437, in gmres
    revcom(b, x, restrt, work, work2, iter_, resid, info, ndx1, ndx2, ijob)
SystemError: NULL result without error in PyObject_Call



Note that I have tested gmres with real sparse matrices, and it runs fine: Indeed

C = sparse.lil_matrix((10, 10))
C.setdiag(rand(10))
C = C.tocsr()
c = rand(10)

x = spla.spsolve(C, c)
print "spsolve: norm(C*x - c) = ", norm(C*x - c)

(x,info) = spla.gmres(C, c)
print "gmres: norm(C*x - c) = ", norm(C*x - c)

gives

spsolve: norm(C*x - c) =  6.93889390391e-18
gmres: norm(C*x - c) =  5.28860261481e-16


I have looked on the web for solutions but haven't found any. Some very old posts indicate similar errors but they don't come with an answer,
and I imagine that if those were due to bugs, they should have been fixed by now...

Am I doing something stupid here, or is that a real problem ? Is somebody aware of a solution ?
(I am using scipy version 0.10.1)

Thanks,
Martin

-- 
Martin Campos Pinto
LJLL, UPMC & CNRS