[SciPy-User] help with scipy cholesky_banded

[Phil Cummins]

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>