Wrapping geqp3 (Rank revealing QR)
Hi all, I was trying to wrap geqp3...(my homework as suggested by Arnd ;-) ) It turns out that it is not straightforward (at least for me). I have used an example taken from generic_flapack.pyf subroutine <tchar=s,d,c,z>geqrf(m,n,a,tau,work,lwork,info) ! qr_a,tau,work,info = geqrf(a,lwork=3*n,overwrite_a=0) ! Compute a QR factorization of a real M-by-N matrix A: ! A = Q * R. callstatement (*f2py_func)(&m,&n,a,&m,tau,work,&lwork,&info) callprotoargument int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int* integer intent(hide),depend(a):: m = shape(a,0) integer intent(hide),depend(a):: n = shape(a,1) <type_in> dimension(m,n),intent(in,out,copy,out=qr) :: a <type_in> dimension(MIN(m,n)),intent(out) :: tau integer optional,intent(in),depend(n),check(lwork>=n||lwork==-1) :: lwork=3*n <type_in> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work integer intent(out) :: info end subroutine <tchar=s,d,c,z>geqrf Some observations: The number of arguments differ wr.t. to complex (C,Z) and real matrices (S,D). How is this handled with scipy ? SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, $ INFO ) SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO ) Here is my first (incomplete !!) attempt to adapt the wrapper . How can I proceed ? ! -*- f90 -*- ! Note: the context of this file is case sensitive. python module wrap_lap ! in interface ! in :wrap_lap subroutine <tchar=s,d>geqp3(m,n,a,lda,jpvt,tau,work,lwork,info) ! in :wrap_lap:dgeqp3.f ! rrqr_a,tau,work,info = geqp3(a,lwork=...,overwrite_a=0) ! Compute the a QR factorization with column pivoting of a matrix A ! A*P = Q*R using Level 3 BLAS callstatement (*f2py_func)(&m,&n,a,&m,tau,work,&lwork,&info) callprotoargument int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int* integer intent(hide),depend(a):: m=shape(a,0) integer intent(hide),depend(a):: n=shape(a,1) <type_in> dimension(m,n), intent(in,out,copy,out=rrqr) :: a <type_in> dimension(MIN(m,n)),intent(out) :: tau integer dimension(*) :: jpvt double precision dimension(*), intent(out) :: tau double precision dimension(*), intent(out) :: work <type_in> dimension(lwork),intent(hide,cache),depend(lwork):: work integer intent(out):: info end subroutine <tchar=s,d>geqp3 end interface end python module wrap_lap ! This file was auto-generated with f2py (version:2_2163). ! See http://cens.ioc.ee/projects/f2py2e/ ~ Nils
On Fri, 24 Feb 2006, Nils Wagner wrote:
Hi all,
I was trying to wrap geqp3...(my homework as suggested by Arnd ;-) ) It turns out that it is not straightforward (at least for me).
I have used an example taken from generic_flapack.pyf
First, all addition to lapack wrapping should be commited to scipy/Lib/lib/lapack. Signatures for geqp3 should be added to flapack_llsc.pyf.src file in scipy/Lib/lib/lapack.
Some observations: The number of arguments differ wr.t. to complex (C,Z) and real matrices (S,D). How is this handled with scipy ?
SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, $ INFO ) SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
See the first signature in flapack_esv.pyf.src that demonstrates how to deal with such a situation. Basically, you should use <_1=,,rwork\,,\2> instead of rwork, in the argument list of <prefix>geqp3. Check out also the corresponing bits in callstatement, callprotoargument are handled in <prefix><sym>ev. HTH, Pearu
Pearu Peterson wrote:
On Fri, 24 Feb 2006, Nils Wagner wrote:
Hi all,
I was trying to wrap geqp3...(my homework as suggested by Arnd ;-) ) It turns out that it is not straightforward (at least for me).
I have used an example taken from generic_flapack.pyf
First, all addition to lapack wrapping should be commited to scipy/Lib/lib/lapack. Signatures for geqp3 should be added to flapack_llsc.pyf.src file in scipy/Lib/lib/lapack.
Some observations: The number of arguments differ wr.t. to complex (C,Z) and real matrices (S,D). How is this handled with scipy ?
SUBROUTINE CGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK, $ INFO ) SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
See the first signature in flapack_esv.pyf.src that demonstrates how to deal with such a situation. Basically, you should use <_1=,,rwork\,,\2> instead of rwork, in the argument list of <prefix>geqp3. Check out also the corresponing bits in callstatement, callprotoargument are handled in <prefix><sym>ev.
HTH, Pearu
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Pearu, Thank you for your quick reply ! The next issues are 1. How do I handle JPVT ? 2. How do I find the size of LWORK in *geqp3 ? I am willing to learn it :-) but the contents of flapack_esv.pyf.src is somehow cryptic in many places e.g. callstatement if(irange_capi==Py_None);else{irange[0]++;irange[1]++;}(*f2py_func)((compute_v?"V":"N"),(vrange_capi==Py_None?(irange_capi==Py_None?"A":"I"):"V"),(lower?"L":"U"),&n,a,&n,vrange,vrange+1,irange,irange+1,&atol,&m,w,z,&ldz,isuppz,work,&lwork,<_2=,,rwork\,&lrwork\,,\2>iwork,&liwork,&info);if(irange_capi==Py_None);else{irange[0]--;irange[1]--;}if(vrange_capi==Py_None);else {capi_w_tmp-\>dimensions[0]=capi_z_tmp-\>dimensions[1]=m;/*capi_z_tmp-\>strides[0]=m*capi_z_tmp-\>descr-\>elsize;*/} Nils
On Fri, 24 Feb 2006, Nils Wagner wrote:
The next issues are
1. How do I handle JPVT ?
2. How do I find the size of LWORK in *geqp3 ?
Read the headers of ?geqp3.f files where the arguments of related subroutines are explained.
I am willing to learn it :-) but the contents of flapack_esv.pyf.src is somehow cryptic in many places e.g.
callstatement if(irange_capi==Py_None);else{irange[0]++;irange[1]++;}(*f2py_func)((compute_v?"V":"N"),(vrange_capi==Py_None?(irange_capi==Py_None?"A":"I"):"V"),(lower?"L":"U"),&n,a,&n,vrange,vrange+1,irange,irange+1,&atol,&m,w,z,&ldz,isuppz,work,&lwork,<_2=,,rwork\,&lrwork\,,\2>iwork,&liwork,&info);if(irange_capi==Py_None);else{irange[0]--;irange[1]--;}if(vrange_capi==Py_None);else {capi_w_tmp-\>dimensions[0]=capi_z_tmp-\>dimensions[1]=m;/*capi_z_tmp-\>strides[0]=m*capi_z_tmp-\>descr-\>elsize;*/}
That is simple C code in one line and takes care of the offset between C and Fortran array indices: if (irange_capi==Py_None) ; else { irange[0]++; irange[1]++; } (*f2py_func)((compute_v?"V":"N"), (vrange_capi==Py_None?(irange_capi==Py_None?"A":"I"):"V"), (lower?"L":"U"), &n,a,&n,vrange,vrange+1,irange,irange+1,&atol,&m, w,z,&ldz,isuppz,work,&lwork, <_2=,,rwork\,&lrwork\,,\2>iwork,&liwork,&info); if (irange_capi==Py_None) ; else { irange[0]--; irange[1]--; } if (vrange_capi==Py_None) ; else { capi_w_tmp-\>dimensions[0]=capi_z_tmp-\>dimensions[1]=m; /*capi_z_tmp-\>strides[0]=m*capi_z_tmp-\>descr-\>elsize;*/ } See also numpy.distutils.from_template.__doc__. Pearu
Pearu Peterson wrote:
On Fri, 24 Feb 2006, Nils Wagner wrote:
The next issues are
1. How do I handle JPVT ?
2. How do I find the size of LWORK in *geqp3 ?
Read the headers of ?geqp3.f files where the arguments of related subroutines are explained.
Sure. I mean the special handling inside the wrapper. I found several ways for different routines to declare _lwork <_lwork=3*n-1,\0,2*n-1,\2> integer optional,intent(in),depend(n) :: lwork=<_lwork> check(<_lwork>\<=lwork) lwork ! <_lwork=(compute_vl||compute_vr)?4*n:3*n,\0,2*n,\2> integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=<_lwork> check(lwork>=<_lwork>) :: lwork ! <_lwork=3*minmn+MAX(2*minmn\,MAX(maxmn\,nrhs)),\0,2*minmn+MAX(maxmn\,nrhs),\2> integer optional,intent(in),depend(nrhs,minmn,maxmn), & check(lwork>=1) & :: lwork=<_lwork> !check(lwork>=<_lwork>) ! <_lwork=n,\0,\0,\0> integer optional,intent(in),depend(n),check(lwork\>=<_lwork>) :: lwork=<_lwork> <ftype> dimension(lwork),intent(hide,cache),depend(lwork) :: work Which way is suitable here ? Nils
I am willing to learn it :-) but the contents of flapack_esv.pyf.src is somehow cryptic in many places e.g.
callstatement if(irange_capi==Py_None);else{irange[0]++;irange[1]++;}(*f2py_func)((compute_v?"V":"N"),(vrange_capi==Py_None?(irange_capi==Py_None?"A":"I"):"V"),(lower?"L":"U"),&n,a,&n,vrange,vrange+1,irange,irange+1,&atol,&m,w,z,&ldz,isuppz,work,&lwork,<_2=,,rwork\,&lrwork\,,\2>iwork,&liwork,&info);if(irange_capi==Py_None);else{irange[0]--;irange[1]--;}if(vrange_capi==Py_None);else {capi_w_tmp-\>dimensions[0]=capi_z_tmp-\>dimensions[1]=m;/*capi_z_tmp-\>strides[0]=m*capi_z_tmp-\>descr-\>elsize;*/}
That is simple C code in one line and takes care of the offset between C and Fortran array indices:
if (irange_capi==Py_None) ; else { irange[0]++; irange[1]++; } (*f2py_func)((compute_v?"V":"N"), (vrange_capi==Py_None?(irange_capi==Py_None?"A":"I"):"V"), (lower?"L":"U"), &n,a,&n,vrange,vrange+1,irange,irange+1,&atol,&m, w,z,&ldz,isuppz,work,&lwork, <_2=,,rwork\,&lrwork\,,\2>iwork,&liwork,&info); if (irange_capi==Py_None) ; else { irange[0]--; irange[1]--; } if (vrange_capi==Py_None) ; else { capi_w_tmp-\>dimensions[0]=capi_z_tmp-\>dimensions[1]=m; /*capi_z_tmp-\>strides[0]=m*capi_z_tmp-\>descr-\>elsize;*/ }
See also numpy.distutils.from_template.__doc__.
Pearu
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
On Fri, 24 Feb 2006, Nils Wagner wrote:
Read the headers of ?geqp3.f files where the arguments of related subroutines are explained.
Sure. I mean the special handling inside the wrapper. I found several ways for different routines to declare _lwork
<_lwork=3*n-1,\0,2*n-1,\2> integer optional,intent(in),depend(n) :: lwork=<_lwork> check(<_lwork>\<=lwork) lwork
! <_lwork=(compute_vl||compute_vr)?4*n:3*n,\0,2*n,\2> integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=<_lwork> check(lwork>=<_lwork>) :: lwork
! <_lwork=3*minmn+MAX(2*minmn\,MAX(maxmn\,nrhs)),\0,2*minmn+MAX(maxmn\,nrhs),\2> integer optional,intent(in),depend(nrhs,minmn,maxmn), & check(lwork>=1) & :: lwork=<_lwork> !check(lwork>=<_lwork>)
! <_lwork=n,\0,\0,\0> integer optional,intent(in),depend(n),check(lwork\>=<_lwork>) :: lwork=<_lwork> <ftype> dimension(lwork),intent(hide,cache),depend(lwork) :: work
Which way is suitable here ?
Read numpy.distutils.from_template.__doc__ about the syntax of <..> blocks. And read ?geqp3.f to find out what are proper _lwork values for each version of geqp3 subroutine. Pearu
Pearu Peterson wrote:
On Fri, 24 Feb 2006, Nils Wagner wrote:
Read the headers of ?geqp3.f files where the arguments of related subroutines are explained.
Sure. I mean the special handling inside the wrapper. I found several ways for different routines to declare _lwork
<_lwork=3*n-1,\0,2*n-1,\2> integer optional,intent(in),depend(n) :: lwork=<_lwork> check(<_lwork>\<=lwork) lwork
! <_lwork=(compute_vl||compute_vr)?4*n:3*n,\0,2*n,\2> integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=<_lwork> check(lwork>=<_lwork>) :: lwork
! <_lwork=3*minmn+MAX(2*minmn\,MAX(maxmn\,nrhs)),\0,2*minmn+MAX(maxmn\,nrhs),\2> integer optional,intent(in),depend(nrhs,minmn,maxmn), & check(lwork>=1) & :: lwork=<_lwork> !check(lwork>=<_lwork>)
! <_lwork=n,\0,\0,\0> integer optional,intent(in),depend(n),check(lwork\>=<_lwork>) :: lwork=<_lwork> <ftype> dimension(lwork),intent(hide,cache),depend(lwork) :: work
Which way is suitable here ?
Read numpy.distutils.from_template.__doc__ about the syntax of <..> blocks. And read ?geqp3.f to find out what are proper _lwork values for each version of geqp3 subroutine.
Pearu
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
numpy.distutils.from_template.__doc__ Traceback (most recent call last): File "<stdin>", line 1, in ? AttributeError: 'module' object has no attribute 'from_template'
dir (numpy.distutils) ['ScipyTest', '_INSTALLED', '__builtins__', '__config__', '__doc__', '__file__', '__name__', '__path__', '__version__', 'ccompiler', 'exec_command', 'log', 'misc_util', 'test', 'unixccompiler']
Am I missing something ? Nils
On Fri, 24 Feb 2006, Nils Wagner wrote:
Read numpy.distutils.from_template.__doc__ about the syntax of <..> blocks. And read ?geqp3.f to find out what are proper _lwork values for each version of geqp3 subroutine.
numpy.distutils.from_template.__doc__ Traceback (most recent call last): File "<stdin>", line 1, in ? AttributeError: 'module' object has no attribute 'from_template'
dir (numpy.distutils) ['ScipyTest', '_INSTALLED', '__builtins__', '__config__', '__doc__', '__file__', '__name__', '__path__', '__version__', 'ccompiler', 'exec_command', 'log', 'misc_util', 'test', 'unixccompiler']
Am I missing something ?
Yes. Have you tried
import numpy.distutils.from_template as m print m.__doc__
or
help('numpy.distutils.from_template')
or $ pydoc numpy.distutils.from_template ? Pearu
On Fri, 24 Feb 2006 09:19:56 -0600 (CST) Pearu Peterson <pearu@scipy.org> wrote:
On Fri, 24 Feb 2006, Nils Wagner wrote:
Read numpy.distutils.from_template.__doc__ about the syntax of <..> blocks. And read ?geqp3.f to find out what are proper _lwork values for each version of geqp3 subroutine.
numpy.distutils.from_template.__doc__ Traceback (most recent call last): File "<stdin>", line 1, in ? AttributeError: 'module' object has no attribute 'from_template'
dir (numpy.distutils) ['ScipyTest', '_INSTALLED', '__builtins__', '__config__', '__doc__', '__file__', '__name__', '__path__', '__version__', 'ccompiler', 'exec_command', 'log', 'misc_util', 'test', 'unixccompiler']
Am I missing something ?
Yes. Have you tried
import numpy.distutils.from_template as m print m.__doc__
or
help('numpy.distutils.from_template')
or
$ pydoc numpy.distutils.from_template
?
Pearu
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Hi Pearu, Thank you for your hint. I will continue my work on Monday. Thank you again Nils
participants (2)
-
Nils Wagner -
Pearu Peterson