[SciPy-dev] newscipy, band matrix routines

[Nils Wagner]

Arnd Baecker wrote:
>Hi,
>
>a while ago I created f2py wrappers for
>- dsbev
>- dsbevd
>- dgbtrf
>- dgbtrs
>- dgbrfs
>which deal with eigenvalue/vector computation,
>LU decomposition, ... for band matrices.
>
>The present version of this, including examples/tests and
>a little bit of documentation is at
>  http://www.physik.tu-dresden.de/~baecker/python/Band2.zip
>
>I would like to get this into newscipy but we need some assistance:
>- should the wrappers be added to generic_flapack.pyf ?
>
>  (+adding the new functions to setup_linalg.py so that the
>   corresponding single precision wrappers can be skipped)
>
>- Are the only scipy specific tricks the use of
>   <tchar=s,d,c,z> and <type_in>
>  to support the different types?
>- For testing/development purpuses, is there
>  a way to put this into a separate sandbox
>  or is it just easier to use all the tools which
>  are provided by Lib/linalg/?
>- it would be great if an expert (Pearu - are you listening? ;-)
>  has a look at the wrapper (either already now, or
>  when the above steps are done) - see also the questions below.
>
>Any further guidance/comments/suggestions are welcome!
>
>Best, Arnd
>
>
>Some more technical questions
>
>Concerning dsbevx:
> a) abstol parameter:
>    from dsbev.f:
>    "Eigenvalues will be computed most accurately when ABSTOL is
>    set to twice the underflow threshold 2*DLAMCH('S'), not zero."
>
>    Maybe the easiest (best?) is to wrap DLAMCH as well (easy)
>    and let the user provide the value, if necessary?
>    ((Even nicer would be to compute this value internally -
>    is this possible with f2py ? If not, I am not sure if adding this
>    is worth the effort...))
>
> b) Presently I suppress the array q
>         Q       (output) DOUBLE PRECISION array, dimension (LDQ, N)
>                 If JOBZ = 'V', the N-by-N orthogonal matrix used in the
>                                reduction to tridiagonal form.
>                 If JOBZ = 'N', the array Q is not referenced.
>    Should one better return this as well and therefore replace
>      double precision dimension(ldq,ldq),intent(hide),depend(ldq) :: q
>    by
>      double precision dimension(ldq,ldq),intent(out),depend(ldq) :: q
>    ?
>
> c) Should one make vl,vu,il,iu, optional as well ?
>    Background: dsbevx allows to compute
>      range: 0   = 'A': all eigenvalues will be found;
>             1   = 'V': all eigenvalues in the half-open interval (vl,vu]
>                        will be found;
>             2   = 'I': the il-th through iu-th eigenvalues will be found.
>    so we may have that either vl,vu or il,iu or neither of
>    the two pairs are used.
>
>    If one makes all vl,vu,il,iu, optional, one could set
>    set the defaults to
>          il=0 iu=10
>          vl=0.0 vu=20.0         # this of course absolutely
>                                 # dependent on the problem,
>                                 # i.e in general nonsense...
>
>
>_______________________________________________
>Scipy-dev mailing list
>Scipy-dev at scipy.net
>http://www.scipy.net/mailman/listinfo/scipy-dev
>  


There is also some progress in improving LAPACK in this context.
http://www.cs.berkeley.edu/~demmel/Sca-LAPACK-Proposal.pdf

Nils