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

Robert Cimrman cimrman3 at ntc.zcu.cz
Thu Dec 14 04:17:32 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)

1. If you use umfpack, k.dtype.char should equal 'D', not 'F' (Look at 
Nils response)
2. also, in case of umfpack usage, use

k=n.empty( (3,3), dtype = 'D' )
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;
k=s.sparse.csr_matrix(k,dtype='D')

because assigning to a CSR matrix via [] (__setitem__) is currently 
buggy. You could use e.g. the following form of the CSR matrix constructor
csr_matrix((data, ij), [dims=(M, N), nzmax=nzmax])
for real data, then, see the doc.

There is now discussion on scipy-dev considering an overhaul of the 
sparse module (speed-wise and functionality-wise and antibug-wise), so 
hopefully bright future awaits us :).

r.


More information about the SciPy-user mailing list