[SciPy-dev] band matrices in scipy

[Arnd Baecker]

Hi,

I am in the progress of wrapping the Lapack dsbev.f routine
(for band-matrices). Thanks to the amazing f2py I got pretty far
(actually in addition the good docs and Fernandos notes helped a lot!).
Once this is finished I would be happy to see this included in scipy ;-).

But before this could happen I have a couple of questions
(sorry, the text below is long,  but I wanted to be explicit ;-):
  - Is the following way of calling dsbev from python alright,
    (eg., I omitted ldab on the python side.)
    Would that parameter be useful to someone so that it is
    better to expose it (the full details of my few pyf lines are at
    the end of this mail)?

    In [7]: MODULENAME.dsbev?
    Type:           fortran
    String Form:    <fortran object at 0x4028ebf0>
    Namespace:      Interactive
    Docstring:
        dsbev - Function signature:
          w,z,info = dsbev(ab,[compute_v,lower])
        Required arguments:
          ab : input rank-2 array('d') with bounds (ldab,*)
        Optional arguments:
          compute_v := 1 input int
          lower := 0 input int
        Return objects:
          w : rank-1 array('d') with bounds (n)
          z : rank-2 array('d') with bounds (ldz,ldz)
          info : int

  - How can one supply better doc-strings to this routine?
    I.e., help(MODULENAME.dsbev) just gives
      "Help on fortran: <fortran object>"
    and scipy.info gives the same as ? from Ipython.

    I would like to add a more detailed description and an example.
    Is this somehow possible to do from within the .pyf wrapper?
        (I know that this is a pure f2py question, so maybe I should
         ask this question as well on the f2py mailing list...)

  - And one more technical aspect (before I forget about it):
    Is the wrapping done such that for _large_ matrices
    no unnecessary copies of arrays are made?

  - One thing which puzzled me a bit (presumably it is
    alright) concerns the return object z which contains the
    eigenvectors if compute_v=1. For compute_v=0
    I changed ldz in the wrapper such that it is of length 1
    (it is never referenced by dsbev.f). So if someone
    is just interested in eigenvalues and not eigenvectors
    a dummy variable for z has to be given.
    Is this alright (for a more "low" level routine like this)?
    In principle one could write a wrapper a la scipy.linalg.eig
    for the band matrices, but for my purposes I don't mind
    about the (sometimes needed) dummy variable.

  - Is the following correct?
    To get this into scipy.linalg.flapack the
    lines from MODULENAME.pyf have to be added to generic_flapack.pyf
    and <tchar=s,d> and <type_in_c>* and <type_in>
    have to be added at the corresponding place. And that's all?

  - Unittests: I think it would be good to have a few unit-tests
    for these routines. Presently I just created a sample
    10x10 band matrix (see below), converted this to upper band storage
    and compared the eigenvalues/eigenvectors  of scipy.linalg.eig
    with those of dsbev.
    In principle one should do this for both upper and lower
    triangle storage and for single and double.

    Would something like that be a reasonable unit-test for inclusion in
scipy
    or is there some better approach
    (better test-matrices etc.)?

  - Further routines I am thinking of wrapping in this context
    a) for symmetric band matrices: dsbevd, dsbevx
    b) for Hermitian band matrices: zhbev, zhbevd, zhbevx
    and their single precision partners.

    Any opinions?

Many thanks,


Arnd

#############################
And here the full example plus test:

MODULENAME.pyf:

# ------ snip here -------------------
!    -*- f90 -*-
python module MODULENAME ! in
    interface  ! in :MODULENAME
        !subroutine dsbev(ab,compute_v,lower) ! in :MODULENAME:dsbev.f
        subroutine dsbev(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,info) !
in :MODULENAME:dsbev.f


          callstatement
(*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,&kd,ab,&ldab,w,z,&ldz,work,&info)

          callprotoargument
char*,char*,int*,int*,double*,int*,double*,double*,int*,double*,int*

          double precision dimension(ldab,*) :: ab

          integer optional,intent(in):: compute_v = 1
          check(compute_v==1||compute_v==0) compute_v
          integer optional,intent(in),check(lower==0||lower==1) :: lower =
0

          ! FIXME:  does one really need the ldab as optional argument ?
          !integer optional,check(shape(ab,0)==ldab),depend(ab) ::
ldab=shape(ab,0)
          integer intent(hide),depend(ab) :: ldab=shape(ab,0)
          integer intent(hide),depend(ab) :: n=shape(ab,1)
          integer intent(hide),depend(ab) :: kd=shape(ab,0)-1

          double precision dimension(n),intent(out),depend(n) :: w
          double  precision dimension(ldz,ldz),intent(out),depend(ldz) ::
z

          ! For compute_v=1 z is used and contains the eigenvectors
          integer intent(hide),depend(n) :: ldz=compute_v?n:1
          double precision dimension(ldz,ldz),depend(ldz) :: z

          double precision dimension(MAX(1,3*n-1)),intent(hide),depend(n)
:: work
          integer intent(out)::info
        end subroutine dsbev
    end interface
end python module MODULENAME

! This file was auto-generated with f2py (version:2.37.233-1545).
! See http://cens.ioc.ee/projects/f2py2e/
! and then modified ...
# ------ snip here -------------------


Then  (dsbev.f from lapack is needed for this)
  f2py --debug-capi -c MODULENAME.pyf dsbev.f -llapack
and a small test is
########### snip here ###########
import MODULENAME
from Numeric import *
from scipy import diag,linalg

N=10
fullmat=diag(-1.0*ones(N-1),-1)+diag(-1.0*ones(N-1),1)+diag(1.0*ones(N))

KD=2        # number of diagonals above the diagonal
LDAB=KD+1
mat=zeros( (LDAB,N),typecode=Float)

for i in arange(LDAB):       # Transform to band matrix
    mat[LDAB-i-1,i:N]=diag(fullmat,i)

w,evec,info=MODULENAME.dsbev(mat)
w_lin,evec_lin=linalg.eig(fullmat)

print sort(w)-sort(w_lin.real)

print evec[:,0]
print evec_lin[:,0]
########### snip here ###########