[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