1

Hi I'm using f2py to wrap the lapack routine dgesvd, by compiling the dgesvd.f file and linking it against llapack, as explained here

according to the docstring, the dgesvd module has the signature:

dgesvd - Function signature:
  dgesvd(jobu,jobvt,m,n,a,s,u,vt,work,lwork,info,[lda,ldu,ldvt])
Required arguments:
  jobu : input string(len=1)
  jobvt : input string(len=1)
  m : input int
  n : input int
  a : input rank-2 array('d') with bounds (lda,*)
  s : input rank-1 array('d') with bounds (*)
  u : input rank-2 array('d') with bounds (ldu,*)
  vt : input rank-2 array('d') with bounds (ldvt,*)
  work : input rank-1 array('d') with bounds (*)
  lwork : input int
  info : input int
Optional arguments:
  lda := shape(a,0) input int
  ldu := shape(u,0) input int
  ldvt := shape(vt,0) input int

Then I use the following ocde to call the module:

mat = rand(20,30)
out_u,out_s,out_vh = zeros((20,20)), zeros((20,)), zeros((30,30))
rows, cols = shape(mat)
workspace = zeros((rows*cols))
out_info = 0

dgesvd(jobu='S', 
    jobvt='S',
    m=rows,
    n=cols,
    a=mat,
    s=out_s,
    u=out_u,
    vt=out_vh,
    work=workspace,
    lwork=rows*cols,
    info=out_info)

Which gives me the correct singular values stored in out_s, but the matrices out_u and out_vh are still only filled with zeros, do I have to do something differently to get the left/right singular vectors as well ?

The code runs through, without any error, which means that out_info is 0.

(The argument 'S' for jobu and jobvt tells the routine to compute only the first min(m,n) singular vectors. Changing it to 'A', does not make any difference)

Any Ideas are highly appreciated! Thanks Mischa

Community
  • 1
  • 1
Mischa Obrecht
  • 2,737
  • 6
  • 21
  • 31
  • Are wrapping the LAPACK routine directly, or writing your own routine that calls the LAPACK routine? – SethMMorton May 27 '12 at 18:15
  • I'm wrapping it directly by using `f2py -L/usr/lib/ -llapack -lblas -dgesvd -c dgesvd.f` The dgesvd routine above is dgesvd.dgesvd. I'm linking against an older lapack version, however. (Source 3.4.1, liblapack: 3.3.1) Hopefully that's not the problem, because I plan to use it on a cluster without recompiling the whole lapack library on it ... – Mischa Obrecht May 27 '12 at 19:58

1 Answers1

5

f2py creates python wrappers to Fortran code, but the python functions that are created are not intended to be called exactly like Fortran code. In Fortran, it is common practice to pass the output variables as an argument to the subroutine. This is not "pythonic"; besides, python doesn't really support subroutines in the same way as Fortran. For this reason, f2py turns your Fortran subroutine into a python function, and thus all output variables are returned by the function, not included in the call signature. So, you would have to call the function this way:

out_s, out_u, out_vh, info = dgesvd(jobu='S', 
                                    jobvt='S',
                                    m=rows,
                                    n=cols,
                                    a=mat,
                                    work=workspace,
                                    lwork=rows*cols)

However, the LAPACK routine is written in FORTRAN77, so it does not have any INTENT declarations for the input/output variables. f2py uses the INTENT declarations to figure out which variables are used as input, and which are to be returned as output. Based on the function signature that you posted, f2py has assumed that all variables are input, which is not what you want. For this reason, I recommend writing your own Fortran 90 wrapper routine that calls dgesvd, so that you can add INTENT declarations yourself to give f2py some hints. I personally would also use the wrapper to allocate the work array to pass to dgesvd so that you don't have to pass it in from python. Exactly how f2py determines the input/output signature is explained here (there are three ways to do it, I prefer the third).

SethMMorton
  • 45,752
  • 12
  • 65
  • 86