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 ###########
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
On Sun, Nov 16, 2003 at 03:50:39PM -0600, Pearu Peterson wrote: [snip]
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.
Looking through pydoc.py (which the builtin help() calls), the problem is that pydoc.py does not recognize a fortranobject as a function-alike. For the interested, look at the method Doc.document. It uses the inspect module, and that only seems to recognize the builtin function-like types. I'll file a bug report asking for something a little smarter. In the meantime, SciPy's info() and IPython's ? magic work just fine. -- Robert Kern kern@ugcs.caltech.edu "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
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
participants (3)
-
Arnd Baecker -
Pearu Peterson -
Robert Kern