1

嗨,我正在使用 f2py 来包装 lapack 例程 dgesvd,方法是编译 dgesvd.f 文件并将其链接到 llapack,如此处所述

根据文档字符串, dgesvd 模块具有签名:

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

然后我使用以下 ocde 调用模块:

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)

这给了我存储在 中的正确奇异值out_s,但是矩阵out_uout_vh仍然只填充零,我是否必须做一些不同的事情来获得左/右奇异向量?

代码运行通过,没有任何错误,out_info即为0。

(jobu 和 jobvt 的参数“S”告诉例程只计算前 min(m,n) 个奇异向量。将其更改为“A”,没有任何区别)

任何想法都受到高度赞赏!谢谢米莎

4

1 回答 1

5

f2py为 Fortran 代码创建 python 包装器,但创建的 python 函数并不打算像 Fortran 代码一样被调用。在 Fortran 中,通常的做法是将输出变量作为参数传递给子例程。这不是“pythonic”;此外,python 并不像 Fortran 那样真正支持子例程。出于这个原因,f2py将您的 Fortran 子例程转换为 python 函数,因此所有输出变量都由函数返回,而不包括在调用签名中。因此,您必须以这种方式调用该函数:

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

但是,LAPACK 例程是用 FORTRAN77 编写的,因此它没有任何INTENT输入/输出变量的声明。 f2py使用INTENT声明来确定哪些变量用作输入,哪些将作为输出返回。根据您发布的函数签名,f2py假设所有变量都是输入的,这不是您想要的。出于这个原因,我建议编写您自己的调用 的 Fortran 90 包装例程dgesvd,以便您可以INTENT自己添加声明以提供f2py一些提示。我个人也会使用包装器来分配要传递的工作数组,dgesvd这样您就不必从 python 传递它。此处f2py解释了如何确定输入/输出签名的确切方式(有三种方法可以做到,我更喜欢第三种)。

于 2012-05-27T20:24:31.927 回答