[SciPy-user] blas/lapack extension

[H Jansen]

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

end function zgeqp3
-------------- 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/437c5d2a/attachment.vcf>