[SciPy-User] Gmres does not work with Sparse Complex Matrices

Martin Campos Pinto mcp.stras@gmail....
Fri Nov 23 00:14:28 CST 2012


Hi all, I am having some troubles using iterative solvers with sparse 
complex matrices.

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



More information about the SciPy-User mailing list