[SciPy-dev] band matrices in scipy

[Arnd Baecker]

Hi,

I continued working on the wrapper for  dsbev, dsbevd and dsbevx.
The present status can be obtained from
  http://www.physik.tu-dresden.de/~baecker/python/Band.zip
Before I continue with the last changes so that it could
be included in scipy it would be nice if someone (Pearu ?;-)
with more knowledge of f2py has a quick look at it.

Concerning dsbevx I have the following questions:
 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...

And one more remark/question concerning the
copying of arrays when using f2py
- sorry if I re-state things already said before, but
I just want to make sure that I understand things correctly here.

First the option -DF2PY_REPORT_ON_ARRAY_COPY=1 is very helpful
(thanks Pearu for pointing this out).
With this I get the following picture (please correct me if I am wrong!):
If one allocates a matrix in python with
   mat=zeros( (N,M))
this is not Fortran contiguous, so a copy is done.
However, declaring
   mat=transpose(zeros(M,N))
one obtains a Fortran contiguous array and no copy is done.
If this is correct, would you say that
the second variant is a good way or is there a better one?

Many thanks,

Arnd