[SciPy-User] Iterative solver for Sparse Complex Matrices

Martin Campos Pinto mcp.stras@gmail....
Fri Nov 23 22:26:31 CST 2012


Hi everybody,

I haven't had much success with my first question so let me ask a new one :

does anybody know an iterative solver (in scipy) that works for linear systems of equations with complex coefficients ?

Thank you,
Martin


On 23 nov. 2012, at 07:14, Martin Campos Pinto <mcp.stras@gmail.com> wrote:

> 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