[SciPy-dev] lapack and blas

Pearu Peterson pearu at cens.ioc.ee
Fri Jan 18 18:04:30 CST 2002


Hi!

Here is a sample pyf file that shows how one can wrap both Fortran LAPACK
and ATLAS LAPACK functions with almost the same signatures. Note that
callstatement is the latest feature in f2py and right now only available
in CVS. The callstatement is quite powerful and allows one to design a
wrapper with almost arbitrary Python signature. Here we needed it because
sgesv functions in lapack and atlas lapack have different signatures
(consider the 'info' argument, for example).

Here are the Python signatures of these wrappers:
  A,ipiv,B,info = f_sgesv(A,B,copy_A=1,copy_B=1)
  A,ipiv,B,info = c_sgesv(A,B,copy_A=1,copy_B=1)

However, a bad news is that lapack and atlas lapack functions are not
compatible. For example,
>>> LU,ipiv,X,info=flapack.c_sgesv([[1,2],[3,4]],[1,1])
>>> print LU
[[ 2.   0.5]
 [ 4.   1. ]]
>>> print X
[-0.5  0.5]
>>> LU,ipiv,X,info=flapack.f_sgesv([[1,2],[3,4]],[1,1])
>>> print LU
[[ 3.          4.        ]
 [ 0.33333334  0.66666669]]
>>> print X
[-0.99999994  0.99999994]

Just thought it would be useful you to know some of the latest features in
f2py and issues in wrapping lapack and atlas lapack functions.

Regards,
	Pearu

!%f90 -*- f90 -*-
python module flapack
interface
   subroutine f_sgesv(n,nrhs,A,lda,ipiv,B,ldb,info)
     fortranname sgesv

     integer depend(A),intent(hide):: n = shape(A,0)
     integer depend(B),intent(hide):: nrhs = shape(B,1) 
     real dimension(n,n),check(shape(A,0)==shape(A,1)),
intent(in,out,copy) :: A
     integer intent(hide),depend(n) :: lda=n
     integer dimension(n),depend(n),intent(out) :: ipiv
     real
dimension(n,nrhs),check(shape(A,0)==shape(B,0)),depend(n),intent(in,out,copy) :: B
     integer intent(hide),depend(n) :: ldb=n
     integer intent(out):: info

   end subroutine sgesv   

   subroutine c_sgesv(n,nrhs,A,lda,ipiv,B,ldb,info)

!
!    Here starts ATLAS specific stuff
!
     fortranname  clapack_sgesv
     intent(c)
     intent(c) :: c_sgesv
     callstatement {extern int clapack_sgesv();  info =clapack_sgesv(102,n,nrhs,A,lda,ipiv,B,ldb); }     
!
! Here ends ATLAS specific stuff, what follows is identical to above
!
     integer depend(A),intent(hide):: n = shape(A,0)
     integer depend(B),intent(hide):: nrhs = shape(B,1) 
     real dimension(n,n),check(shape(A,0)==shape(A,1)),
intent(in,out,copy) :: A
     integer intent(hide),depend(n) :: lda=n
     integer dimension(n),depend(n),intent(out) :: ipiv
     real
dimension(n,nrhs),check(shape(A,0)==shape(B,0)),depend(n),intent(in,out,copy) :: B
     integer intent(hide),depend(n) :: ldb=n
     integer intent(out):: info

   end subroutine c_sgesv
end interface
end python module flapack




More information about the Scipy-dev mailing list