[SciPy-user] wrapping lapack
H Jansen
h.jansen at fel.tno.nl
Tue Jun 3 12:49:13 EDT 2003
<< earlier message was posted in the wrong thread >>
I would appreciate some help on the following:
In need for a pivoting QR decomposition of a rectangular matrix, I set
out to extend scipy.linalg.lapack by writing f2py the wrapper code for
the <tchar=s,d,c,z>geqp3(...). Having mastered the basic f2py skills of
writing python extensions for C code (I'm using ATLAS with full lapack
library), I discovered that the (optimized) ATlAS/lapack routines are
prefixed with "clapack_" (e.g. clapack_gesv), which in turn are wrappers
for "ATL_" (ATLAS) optimized C functions. According to this philosophy,
for each lapack function two more layers are therefore needed between a
python wrapper function and the original lapack C code: one
"clapack_"-like wrapper and one "ATL_"-like wrapper. Am I correct? Is
this what the scipy authors have in mind when planning further
extensions of scipy.linalg?
Or can it be done simpler ...? How about the storage order then?
Another question: why are there two "*.pyf" files: a clapack.pyf file
and a generic_clapack.pyf file?
I've appended my _geqp3 code for the clapack.pyf and generic_clapack.pyf
files.
>>>>>>>>>>>>>>>>>>>>> generic_clapack.pyf <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
function <tchar=s,d,c,z>geqp3(m,n,a,piv,tau,work,lwork,info,rowmajor)
! qr,tau,piv,info = geqp3(a, piv=0, rowmajor=1, overwrite_a=0)
!
! Computes a QR factorization with column pivoting of matrix a so
that
! A*P = Q*R
! using Level 3 BLAS.
! A (input/output) REAL array, dimension (LDA,N)
! On entry, the M-by-N matrix A.
! On exit, the upper triangle of the array contains the
! min(M,N)-by-N upper trapezoidal matrix R; the elements below
! the diagonal, together with the array TAU, represent the
! orthogonal matrix Q as a product of min(M,N) elementary
! reflectors.
! piv pivots columns (input/output) INTEGER array, dimension (N)
! On entry:
! if JPVT(J).ne.0, the J-th column of A is permuted to
! the front of A*P (a leading column);
! if JPVT(J)=0, the J-th column of A is a free column.
! On exit:
! if JPVT(J)=K, then the J-th column of A*P was the the
! K-th column of A.
! TAU (output) REAL array, dimension (min(M,N))
! The scalar factors of the elementary reflectors.
fortranname <tchar=s,d,c,z>geqp3_
integer intent(c,hide) :: <tchar=s,d,c,z>geqp3
callstatement <tchar=s,d,c,z>geqp3_return_value = info =
(*f2py_func)
(102-rowmajor, m, n,
a, m, piv, tau, work, lwork)
callprotoargument const int,const int,const int,<type_in_c>*,const
int,int*,<type_in_c>*,<type_in_c>*,const int
integer optional,intent(in),check(rowmajor==1||rowmajor==0) ::
rowmajor = 1
integer depend(a),intent(hide):: m = shape(a,0)
integer depend(a),intent(hide):: n = shape(a,1)
<type_in> dimension(m,n) :: a
integer dimension(n),depend(n),intent(in,out) :: piv
integer depend(work),intent(hide),check(lwork>=3*n+1):: lwork =
shape(work,0)
<type_in> dimension(lwork),intent(in) :: work
integer intent(out)::info
<type_in> dimension(MIN(m,n)),intent(out,hide) :: tau
intent(c,in,out,copy,out=qr) a
end function <tchar=s,d,c,z>geqp3
>>>>>>>>>>>>>>>>>>>>> clapack.pyf (s,d,c and z codes) <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
function sgeqp3(m,n,a,piv,tau,work,lwork,info,rowmajor)
! qr,tau,piv,info = geqp3(a, piv=0, rowmajor=1, overwrite_a=0)
!
! Computes a QR factorization with column pivoting of matrix a so
that
! A*P = Q*R
! using Level 3 BLAS.
! A (input/output) REAL array, dimension (LDA,N)
! On entry, the M-by-N matrix A.
! On exit, the upper triangle of the array contains the
! min(M,N)-by-N upper trapezoidal matrix R; the elements below
! the diagonal, together with the array TAU, represent the
! orthogonal matrix Q as a product of min(M,N) elementary
! reflectors.
! piv pivots columns (input/output) INTEGER array, dimension (N)
! On entry:
! if JPVT(J).ne.0, the J-th column of A is permuted to
! the front of A*P (a leading column);
! if JPVT(J)=0, the J-th column of A is a free column.
! On exit:
! if JPVT(J)=K, then the J-th column of A*P was the the
! K-th column of A.
! TAU (output) REAL array, dimension (min(M,N))
! The scalar factors of the elementary reflectors.
fortranname sgeqp3
integer intent(c,hide) :: sgeqp3
callstatement sgeqp3_return_value = info =
(*f2py_func)(102-rowmajor,m,n,a,m,piv,tau,work,lwork)
callprotoargument const int,const int,const int,<type_in_c>*,const
int,int*,<type_in_c>*,<type_in_c>*,const int
integer optional,intent(in),check(rowmajor==1||rowmajor==0) ::
rowmajor = 1
integer depend(a),intent(hide):: m = shape(a,0)
integer depend(a),intent(hide):: n = shape(a,1)
<type_in> dimension(m,n) :: a
integer dimension(n),depend(n),intent(in,out) :: piv
integer depend(work),intent(hide),check(lwork>=3*n+1):: lwork =
shape(work,0)
<type_in> dimension(lwork),intent(in) :: work
integer intent(out)::info
integer dimension(MIN(m,n)), intent(out) :: tau
intent(c,in,out,copy,out=qr) a
end function sgeqp3
function dgeqp3(m,n,a,piv,tau,work,lwork,info,rowmajor)
! qr,tau,piv,info = geqp3(a, piv=0, rowmajor=1, overwrite_a=0)
!
! Computes a QR factorization with column pivoting of matrix a so
that
! A*P = Q*R
! using Level 3 BLAS.
! A (input/output) REAL array, dimension (LDA,N)
! On entry, the M-by-N matrix A.
! On exit, the upper triangle of the array contains the
! min(M,N)-by-N upper trapezoidal matrix R; the elements below
! the diagonal, together with the array TAU, represent the
! orthogonal matrix Q as a product of min(M,N) elementary
! reflectors.
! piv pivots columns (input/output) INTEGER array, dimension (N)
! On entry:
! if JPVT(J).ne.0, the J-th column of A is permuted to
! the front of A*P (a leading column);
! if JPVT(J)=0, the J-th column of A is a free column.
! On exit:
! if JPVT(J)=K, then the J-th column of A*P was the the
! K-th column of A.
! TAU (output) REAL array, dimension (min(M,N))
! The scalar factors of the elementary reflectors.
fortranname dgeqp3
integer intent(c,hide) :: dgeqp3
callstatement dgeqp3_return_value = info =
(*f2py_func)(102-rowmajor,m,n,a,m,piv,tau,work,lwork)
callprotoargument const int,const int,const int,<type_in_c>*,const
int,int*,<type_in_c>*,<type_in_c>*,const int
integer optional,intent(in),check(rowmajor==1||rowmajor==0) ::
rowmajor = 1
integer depend(a),intent(hide):: m = shape(a,0)
integer depend(a),intent(hide):: n = shape(a,1)
<type_in> dimension(m,n) :: a
integer dimension(n),depend(n),intent(in,out) :: piv
integer depend(work),intent(hide),check(lwork>=3*n+1):: lwork =
shape(work,0)
<type_in> dimension(lwork),intent(in) :: work
integer intent(out)::info
integer dimension(MIN(m,n)),intent(out,hide) :: tau
intent(c,in,out,copy,out=qr) a
end function dgeqp3
function cgeqp3(m,n,a,piv,tau,work,lwork,info,rowmajor)
! qr,tau,piv,info = geqp3(a, piv=0, rowmajor=1, overwrite_a=0)
!
! Computes a QR factorization with column pivoting of matrix a so
that
! A*P = Q*R
! using Level 3 BLAS.
! A (input/output) REAL array, dimension (LDA,N)
! On entry, the M-by-N matrix A.
! On exit, the upper triangle of the array contains the
! min(M,N)-by-N upper trapezoidal matrix R; the elements below
! the diagonal, together with the array TAU, represent the
! orthogonal matrix Q as a product of min(M,N) elementary
! reflectors.
! piv pivots columns (input/output) INTEGER array, dimension (N)
! On entry:
! if JPVT(J).ne.0, the J-th column of A is permuted to
! the front of A*P (a leading column);
! if JPVT(J)=0, the J-th column of A is a free column.
! On exit:
! if JPVT(J)=K, then the J-th column of A*P was the the
! K-th column of A.
! TAU (output) REAL array, dimension (min(M,N))
! The scalar factors of the elementary reflectors.
fortranname cgeqp3
integer intent(c,hide) :: cgeqp3
callstatement cgeqp3_return_value = info =
(*f2py_func)(102-rowmajor,m,n,a,m,piv,tau,work,lwork)
callprotoargument const int,const int,const int,<type_in_c>*,const
int,int*,<type_in_c>*,<type_in_c>*,const int
integer optional,intent(in),check(rowmajor==1||rowmajor==0) ::
rowmajor = 1
integer depend(a),intent(hide):: m = shape(a,0)
integer depend(a),intent(hide):: n = shape(a,1)
<type_in> dimension(m,n) :: a
integer dimension(n),depend(n),intent(in,out) :: piv
integer depend(work),intent(hide),check(lwork>=3*n+1):: lwork =
shape(work,0)
<type_in> dimension(lwork),intent(in) :: work
integer intent(out)::info
integer dimension(MIN(m,n)),intent(out,hide) :: tau
intent(c,in,out,copy,out=qr) a
end function cgeqp3
function zgeqp3(m,n,a,piv,tau,work,lwork,info,rowmajor)
! qr,tau,piv,info = geqp3(a, piv=0, rowmajor=1, overwrite_a=0)
!
! Computes a QR factorization with column pivoting of matrix a so
that
! A*P = Q*R d
! using Level 3 BLAS.
! A (input/output) REAL array, dimension (LDA,N)
! On entry, the M-by-N matrix A.
! On exit, the upper triangle of the array contains the
! min(M,N)-by-N upper trapezoidal matrix R; the elements below
! the diagonal, together with the array TAU, represent the
! orthogonal matrix Q as a product of min(M,N) elementary
! reflectors.
! piv pivots columns (input/output) INTEGER array, dimension (N)
! On entry:
! if JPVT(J).ne.0, the J-th column of A is permuted to
! the front of A*P (a leading column);
! if JPVT(J)=0, the J-th column of A is a free column.
! On exit:
! if JPVT(J)=K, then the J-th column of A*P was the the
! K-th column of A.
! TAU (output) REAL array, dimension (min(M,N))
! The scalar factors of the elementary reflectors.
fortranname zgeqp3
integer intent(c,hide) :: zgeqp3
callstatement zgeqp3_return_value = info =
(*f2py_func)(102-rowmajor,m,n,a,m,piv,tau,work,lwork)
callprotoargument const int,const int,const int,<type_in_c>*,const
int,int*,<type_in_c>*,<type_in_c>*,const int
integer optional,intent(in),check(rowmajor==1||rowmajor==0) ::
rowmajor = 1
integer depend(a),intent(hide):: m = shape(a,0)
integer depend(a),intent(hide):: n = shape(a,1)
<type_in> dimension(m,n) :: a
integer dimension(n),depend(n),intent(in,out) :: piv
integer depend(work),intent(hide),check(lwork>=3*n+1):: lwork =
shape(work,0)
<type_in> dimension(lwork),intent(in) :: work
integer intent(out)::info
integer dimension(MIN(m,n)),intent(out,hide) :: tau
intent(c,in,out,copy,out=qr) a
-------------- next part --------------
A non-text attachment was scrubbed...
Name: h.jansen.vcf
Type: text/x-vcard
Size: 479 bytes
Desc: Card for H Jansen
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20030603/1f9043fc/attachment.vcf>
More information about the SciPy-User
mailing list