[SciPy-dev] newscipy, band matrix routines

[Arnd Baecker]

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...