[SciPy-User] help with scipy cholesky_banded
Phil Cummins
phil.cummins@anu.edu...
Mon Aug 23 18:22:29 CDT 2010
Hi,
I am using Enthought epd64 version 6.2-2 on Mac OS 10.6.4. The version of scipy used is 0.8.0.dev6485.
I am trying to create a finite element solution to a simple, 1D differential equation. I have constructed a stiffness matrix and called scipy.linalg.cholesky_banded to decompose it. Cholesky_banded complains that the matrix is not positive definite, which should not happen if it is a properly constructed FEM stiffness matrix. I checked my construction many times, and finally started looking into why cholesky_banded thinks it's not positive definite
My understanding is that, since my matrix is real, cholesky_banded reduces to a call to the lapack routine dpbtrf, and since I don't think I'm using the 'blocked' version, I think it should be using dpbtrf2. So I went and found dpbtrf2, and discovered that its check for positive definiteness appears to consist of testing whether any of the diagonals are <= 0. I downloaded dpbtrf2 and commented out the lines that called external routines (i.e., leaving pretty much just the check for positive definiteness). I used f2py to import it as a module, so that I could investigate why the positive definitiveness check fails.
Attached is a python script which sets the elements of the stiffness matrix. It then does the following (note that the matrix Stiff is banded symmetric, with only 1 super diagonal, stored as upper triangular):
import scipy.linalg
import math
from numpy import *
import dpbtf2 # My own dpbtrf2 module created with "f2py -c dpbtf2.pyf dpbtf2.f"
.
.
.
for i in range(0,len(rg)):
if Stiff[1][i] <= 0.:
print 'non positive diagonal %d <= 0.: %g' % (i,Stiff[1][i])
info = dpbtf2.dpbtf2('U',len(rg),1,Stiff,2)
print 'dpbtf2:info=',info
C, info = scipy.linalg.flapack.dpbtrf(Stiff)
print 'dpbtrf:info=',info
So, in python I check whether any diagonal elements are <= 0., then I check this with my modified dpbtrf2 fortran subroutine, and then I call the linalg one. Here's what happens:
$ ./scipy_bug.py
dpbtf2:info= 0
dpbtrf:info= 1001
Python doesn't find a diagonal element <=0., nor does the fortran routine dpbtrf2. But the scipy.linalg.dpbtrf does - so the latter fails.
Can anyone please help me understand why scipy.linalg.dpbtr doesn't seem to work for me?
Thanks,
- Phil
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: scipy_bug.py
Type: text/x-python-script
Size: 1676 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment.bin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dpbtf2.f
Type: application/octet-stream
Size: 5795 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment.obj
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment-0002.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: dpbtf2.pyf
Type: application/octet-stream
Size: 627 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment-0001.obj
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100823/0bd572bd/attachment-0003.html
More information about the SciPy-User
mailing list