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@scipy.org
http://projects.scipy.org/mailman/listinfo/scipy-user