[SciPy-Dev] Gmres with Sparse Complex Matrices gives Segmentation fault (or NULL result)
Martin Campos Pinto
mcp.stras@gmail....
Fri Nov 23 23:50:51 CST 2012
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
More information about the SciPy-Dev
mailing list