通过f2py编译后使用expokit.zgexpv

Question

我想计算表达式，比如exp（At）v，其中A是复数，scipy.sparse.csr矩阵，t是实数，v是复数的一个numpy一维数组

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

expokit.f中有一个不同函数（expokit.zgpadm）的[test][2]文件，它在我的Ubuntu机器上按预期执行，但我对expokit.zgexpv函数感兴趣

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

>>> 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)