问题描述
我想计算表达式,例如 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
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)