[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