通过 f2py 编译后使用 expokit.zgexpv

Question

我想计算表达式，例如 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,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,t,v,w,tol,nrmA,wsp,iwsp,matvec,itrace,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

Answer 1

expokit 的 matvec 子程序调用如下

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

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

我希望此评论能帮助某人真正找到答案。

---

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

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