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.
_______________________________________________
SciPy-user mailing list
SciPy-user@[...].org
http://projects.scipy.org/mailman/listinfo/scipy-user