Indefinite coefficient matrix
Hi all, I wan't to solve a system of linear equations A x = b, where A is symmetric but indefinite. Thus one cannot use a Cholesky factorization L L^T. AFAIK one can use an L D L^T - factorization, where D is a diagonal matrix. Is this factorization available in numpy/scipy ? Nils
Nils Wagner wrote:
Hi all,
I wan't to solve a system of linear equations
A x = b,
where A is symmetric but indefinite. Thus one cannot use a Cholesky factorization L L^T. AFAIK one can use an L D L^T - factorization, where D is a diagonal matrix. Is this factorization available in numpy/scipy ?
Nils, you can use umfpack to solve this. Thanks to Nathan Bell, you can get the actual L, U factors if you need them, too. r.
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)
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@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
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.
participants (3)
-
Mark Starnes
-
Nils Wagner
-
Robert Cimrman