Looking at the source for scipy.linalg.solveh_banded, it just wraps Lapack pbsv. I'm looking for a more efficient solver for tridiagonal (Hermitian, or in my case real symmetric) systems which I think should be provided by the Lapack ptsv function. In addition, solveh_banded will crash if I have too large a dynamic range of (all positive) values along the main diagonal even though this should not be an actual problem (I'm guessing that roundoff makes the smallest values seem effectively negative, so it's seen as having negative eigenvalues) and there's some chance a tridiagonal-specific routine would not hit this issue.
From my reading about Lapack it seems like ptsv should be included in any distribution that has pbsv (documentation always lists them together). I'm not really certain which would be more efficient (pbsv assumes symmetric but with arbitrary bandwidth, while ptsv assumes tridiagonal but not necessary symmetric) but it seemed worth trying ptsv.
Unfortunately, ptsv doesn't seem to be exposed in scipy, which I believe in practice means that it's not included in scipy.linalg.flapack, and is therefore not available via scipy.linalg.get_lapack_funcs(('ptsv',)).
I realize that the Fortran/Lapack linkage with scipy is complicated, but does anybody know why pbsv and ptsv would be treated differently? Is there something I can hand-edit to try wrapping ptsv like pbsv (unfortunately flapack seems to be provided just as a ".so" so I ran into a dead end)?
FWIW I'm using Enthought EPD with the Intel MKL. However, given that scipy.linalg (independent of distribution) always includes solveh_banded, but does not have a tridiagonal solver I'm thinking it must be deeper than just an EPD/MKL issue.
