通过 f2py 编译后使用 expokit.zgexpv

如何解决通过 f2py 编译后使用 expokit.zgexpv

我想计算表达式，例如 exp(At)v，其中 A 是一个复数 scipy.sparse.csr 矩阵，t 是一个实数，v 是一个 numpy 一维复数数组。

目前我有一个文件 [expokit.f][1]，它在第 2420 行附近有一个子程序 ZGEXPV（见链接）。它还具有一些 f2py 意图内涵。我使用 f2py -c expokit.f -m expokit --link-blas_opt --link-lapack_opt 编译了 expokit.f 并且编译没有出错。

expokit.f 中有一个用于不同函数 (expokit.zgpadm) 的 [test][2] 文件，它在我的 Ubuntu 机器上按预期运行良好，但我对 expokit.zgexpv 函数感兴趣.

当我在 f2py 编译后在 Python 中 import expokit 并执行 numpy.info(expokit.zgexpv) 时，我得到以下输出：

>>> np.info(expokit.zgexpv)
zgexpv(m,t,v,w,tol,anorm,wsp,iwsp,matvec,itrace,iflag,[n,lwsp,liwsp,matvec_extra_args])

Wrapper for ``zgexpv``.

Parameters
----------
m : input int
t : input float
v : input rank-1 array('D') with bounds (n)
w : in/output rank-1 array('D') with bounds (n)
tol : in/output rank-0 array(float,'d')
anorm : input float
wsp : input rank-1 array('D') with bounds (lwsp)
iwsp : input rank-1 array('i') with bounds (liwsp)
matvec : call-back function
itrace : input int
iflag : in/output rank-0 array(int,'i')

Other Parameters
----------------
n : input int,optional
    Default: len(v)
lwsp : input int,optional
    Default: len(wsp)
liwsp : input int,optional
    Default: len(iwsp)
matvec_extra_args : input tuple,optional
    Default: ()

Notes
-----
Call-back functions::

  def matvec(n,e_wsp_j1v_n_err,e_wsp_j1v_err): return
  required arguments:
    n : input int
    e_wsp_j1v_n_err : input complex
    e_wsp_j1v_err : input complex

现在我的目标是能够将输入传递给函数如下（python3）：

n=2000
m=30 #(m denotes the number of Krylov steps in the expokit.zgexpv )
A=scipy.sparse.rand(n,n,density=0.01) + 1.0j  * scipy.sparse.rand(n,density=0.01)
v=np.random.rand(n) + 1.0j * np.random.rand(n) #the vector v on which exp(At) acts.

t=1.0
tol=1e-07
itrace=np.array([0])
iflag=np.array([0])
MXLIM=200000
iwsp=np.zeros(MXLIM).astype(int) #sufficiently large array
wsp=np.zeros(MXLIM) #sufficiently large array
nrmA=scipy.linalg.norm(A)

def matvec(n,x,y):
    y = A @ x

answer=expokit.zgexpv(m,nrmA,iflag) #answer should be exp(At)v

显然我在这里做错了。如果有人可以更正我上面提到的代码，那就太好了。我几乎可以肯定我的 matvec 函数有问题，可能还有其他问题。 [1]：https://drive.google.com/file/d/1tvsbBqhsonkVgUbzX0jd5a5uABk2YgS9/view?usp=sharing [2]：https://drive.google.com/file/d/1hqQjzN_bvT-P-fKFzxwJYcaAgqgX84aj/view?usp=sharing

解决方法

问题出现在expokit的matvec子程序如下调用

call matvec(n,wsp(j1v-n),wsp(j1v) )

因此，f2py 仅向您自己的 (python) matvec 函数提供一个整数、一个复数和另一个复数。注意那些复数是标量而不是长度为 n 的向量。

我希望这篇评论能帮助人们真正找到答案。

您的代码中还有一个小错误。 scipy.linalg.norm 不适用于稀疏矩阵。你需要调用以下

nrmA=scipy.sparse.linalg.norm(A)