[SciPy-dev] band matrices in scipy

[Pearu Peterson]

Hi,

On Sun, 16 Nov 2003, Arnd Baecker wrote:

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

Great!

> 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

I would set ldab=kd+1 and hide it from Python. Rationale:
(i) wrappers to other lapack functions do the same with
arguments of similar type, so, why should dsbev be different.
(ii) my guess is that users who need ldab exposed would not
call dsbev from Python in the first place.

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

Currently one cannot.. f2py does not support user specified documentation 
strings yet. However, it is certainly in my TODO list and some 
relevant issues has been already resolved (I am referring to f2py 
multi-line blocks that are prerequisites for implementing 
user-docs-in-pyf-file feature).

However, help(MODULENAME.dsbev) could certainly give more
information than just '<fortran object>'. The reason why it
does not may be due to the fact that a <fortran object> does not
have a __doc__ attribute, instead, when one tries to access
<fortran object>.__doc__ then getattr function is called
that generates documentation on fly. This is f2py feature to 
reduce the size of extension modules.

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

Yes. I have explained many times what is the underlying
mechanism for passing Python objects to Fortran and how
one should design .pyf files in order to minimize the number of 
array copies and at the same time simplify the usage of
wrappers.
Btw, with the recent f2py one can define F2PY_REPORT_ON_ARRAY_COPY 
CPP macro in order to get statistics on copies that happen inside
f2py generated wrappers.

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

Some other pyf-signatures in scipy do the same thing (i.e. set 
ldz=(compute_v?n:1)). So, it is alright for low-level functions
but for high-level functions  one must implement a Python wrapper to get 
rid of z when compute_v=0.

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

Yep. But see also setup_linalg.py that provides a possibility to
skip single precision wrappers (resulting 2 times smaller extension 
modules). So, information for new functions should be added there too.

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

Yes. It assumes that, say, scipy.linalg.eig gives correct results and 
this is checked (in principle) in the corresponding unit-tests.
See the existing scipy.linalg unittests for examples.

>     or is there some better approach
>     (better test-matrices etc.)?

Regards,
Pearu