[SciPy-user] csc, csr sparse complex matrix problems.

Nils Wagner nwagner at iam.uni-stuttgart.de
Thu Dec 14 02:18:12 CST 2006


Mark Starnes wrote:
> Hi everyone,
>
> I've been trying to get this simple example to work for a couple of days 
> now, to no avail.  I noticed at 
> http://projects.scipy.org/scipy/scipy/ticket/329 that there seem to be 
> come problems with csc and complex values but I'm hoping they're 
> unrelated and someone can help me sort this problem!  I couldn't get the 
> csc matrix to work at all here - the csr matrix was having a go and 
> getting the right answer but now fails too!  Anyway, here's the code.  
> I'm using Python2.4.3 on winxp with Scipy0.5.2, Numpy 1.0.1.
>
> Any help will be appreciated.
>
> Best regards,
>
> Mark Starnes.
>
> ***
>
> import numpy as n
> import scipy as s
> from scipy.linsolve import spsolve
> from scipy import linalg as la
> from numpy import zeros,complex64
>
> # solve kx=f for x, where
> #
> #/                    \  /    \    /       \
> #| 1+1i  2+2i  3+3i   | |   2  |   | -5+3i  |
> #|                    | |      |   |        |
> #| 2     3     4+4i   | | 3+2i | = | 1-6i   |
> #|                    | |      |   |        |
> #| 4+9i 5+10i  2+11i  | |  -3  |   | -3+25i |
> #\                    /  \    /    \       /
> #
>
>
> k0=n.array([[1+1j,2+2j,3+3j],[2,3,4+4j],[4+9j,5+10j,2+11j]])
> f0=[-5+3j,1-6j,-3+25j]
>
> # Test numpy routine.
>
> k=n.mat(k0)
> f=n.mat(f0).T
> fsol=s.dot(la.inv(k),f)
> print fsol
> print k*fsol
> # result should = [2, 3+2j, -3].  Numpy passes.
>
> # Test Scipy sparse routines.
> k=s.sparse.csc_matrix((3,3),dtype='F')
> k[0,0]=1+1j ; k[0,1]=2+2j  ; k[0,2]=3+3j ;
> k[1,0]=2    ; k[1,1]=3     ; k[1,2]=4+4j ;
> k[2,0]=4+9j ; k[2,1]=5+10j ; k[2,2]=2+11j;
>
> f=zeros((3),complex64)
> f[0]=-5+3j
> f[1]=1-6j
> f[2]=-3+25j
>
> #fsol2=spsolve(k,f)   # this crashes python with typeD, AND with type F
> #print fsol2
> #print k.dot(fsol2)
>
> k=s.sparse.csr_matrix((3,3),dtype='F')
> k[0,0]=1+1j ; k[0,1]=2+2j  ; k[0,2]=3+3j ;
> k[1,0]=2    ; k[1,1]=3     ; k[1,2]=4+4j ;
> k[2,0]=4+9j ; k[2,1]=5+10j ; k[2,2]=2+11j;
>
> fsol3=spsolve(k,f)   # this crashes python with typeD, originally ok 
> with type F
> print fsol3         
> print k.dot(fsol3)  
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>   
Hi Mark,

This seems to be a an UMFPACK problem. It works fine for me with SuperLU.
I am using the latest svn versions and UMFPACK v4.4

Nils

With UMFPACK

python -i starnes.py
x [[ 2. +1.11022302e-16j]
 [ 3. +2.00000000e+00j]
 [-3. -7.28583860e-17j]]
f [[-5. +3.j]
 [ 1. -6.j]
 [-3.+25.j]]
Traceback (most recent call last):
  File "starnes.py", line 47, in ?
    fsol2=spsolve(k,f)   #
  File "/usr/lib64/python2.4/site-packages/scipy/linsolve/linsolve.py",
line 74, in spsolve
    return umf.linsolve( umfpack.UMFPACK_A, mat, b, autoTranspose = True )
  File
"/usr/lib64/python2.4/site-packages/scipy/linsolve/umfpack/umfpack.py",
line 561, in linsolve
    self.numeric( mtx )
  File
"/usr/lib64/python2.4/site-packages/scipy/linsolve/umfpack/umfpack.py",
line 374, in numeric
    self.symbolic( mtx )
  File
"/usr/lib64/python2.4/site-packages/scipy/linsolve/umfpack/umfpack.py",
line 356, in symbolic
    raise RuntimeError, '%s failed with %s' % (self.funs.symbolic,
RuntimeError: <function umfpack_zi_symbolic at 0x2aaaadff7500> failed
with UMFPACK_ERROR_invalid_matrix

Without UMFPACK

python -i starnes.py
x [[ 2. +1.11022302e-16j]
 [ 3. +2.00000000e+00j]
 [-3. -7.28583860e-17j]]
f [[-5. +3.j]
 [ 1. -6.j]
 [-3.+25.j]]
Use minimum degree ordering on A'+A.
x [ 2. -7.58729969e-17j  3. +2.00000000e+00j -3. -2.77108197e-16j]
f [-5. +3.j  1. -6.j -3.+25.j]
Use minimum degree ordering on A'+A.
x [ 2. -3.44633288e-16j  3. +2.00000000e+00j -3. -2.91704595e-16j]
f [-5. +3.j  1. -6.j -3.+25.j]


Nils
 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: starnes.py
Type: text/x-python
Size: 1436 bytes
Desc: not available
Url : http://projects.scipy.org/pipermail/scipy-user/attachments/20061214/4c6eebdc/attachment.py 


More information about the SciPy-user mailing list