从 Python 调用时，Fortran 代码提供错误的结果

Question

为了提高金融问题的执行速度，我在 Fortran 中编写了核心数字部分，只在 Python 中进行文件访问等。我使用 f2py 编译，并调用子程序fit（见我的帖子下方）

vec=num.num.fit(pfmat,lqmat,20,0.01,20,0.01,pfmat.shape[0],lqmat.shape[0])

我已经测试了 Fortran 代码，其中没有错误，但是，在使用 f2py 编译并从 Python 调用它之后，输出无法重现。有时，运行会返回正确的结果（与 Python 中编程的相同算法相比），其他时候调用返回带有 nan 的向量，有时结果大约为 10^300。目前我不知道如何确定根本原因，非常感谢任何帮助！

subroutine fit(pf,m,lq,n,ns,ds,nv,dv,alloc)
integer, intent(in) :: m
integer, intent(in) :: n
double precision, intent(in) :: pf(m,14)
double precision, intent(in) :: lq(n,14)
integer, intent(in) :: ns
double precision, intent(in) :: ds
integer, intent(in) :: nv
double precision, intent(in) :: dv
double precision :: vit(1,(2*ns+1)*(2*nv+1))
double precision :: ml((2*ns+1)*(2*nv+1),n+1)
double precision, intent(out) :: alloc(1,n+1)
double precision :: matinv(n+1,n+1)
double precision :: mat_post((2*ns+1)*(2*nv+1),n+1)
integer :: i,j,k !
double precision :: f, sig, lm, strike, d1, d2, v0

do i=1,m
  strike=pf(i,1)
  f=pf(i,3)
  lm=log(f/strike)
  sig=pf(i,8)+pf(i,9)*lm+pf(i,10)*lm**2+pf(i,11)*lm**3+pf(i,12)*lm**4+ &
      pf(i,13)*lm**5+pf(i,14)*abs(lm)
  d1=(lm+0.5*sig**2*pf(i,4))/(sig*sqrt(pf(i,4)))
  d2=d1-sig*sqrt(pf(i,4))
  if (pf(i,2)==0)then
    v0=pf(i,5)*pf(i,6)*pf(i,7)*(f*cnd(d1)-strike*cnd(d2))
  else
    v0=pf(i,5)*pf(i,6)*pf(i,7)*(strike*cnd(-d2)-f*cnd(-d1))
  end if
  ! loop over all gridpoints for (u/l shock and vol shock)
  do j=-ns,ns
    do k=-nv,nv
        f=pf(i,3)*(1+j*ds)
        lm=log(f/strike)
        sig=pf(i,8)+pf(i,9)*lm+pf(i,10)*lm**2+pf(i,11)*lm**3+pf(i,12)*lm**4+ &
            pf(i,13)*lm**5+pf(i,14)*abs(lm)
        sig=sig+pf(i,8)*k*dv
        d1=(lm+0.5*sig**2*pf(i,4))/(sig*sqrt(pf(i,4)))
        d2=d1-sig*sqrt(pf(i,4))
        if (pf(i,2)==0)then
          vit(1,(j+ns)*(2*nv+1)+k+nv+1)=vit(1,(j+ns)*(2*nv+1)+k+nv+1)+ &
            pf(i,5)*pf(i,6)*pf(i,7)*(f*cnd(d1)-strike*cnd(d2))-v0
        else
          vit(1,(j+ns)*(2*nv+1)+k+nv+1)=vit(1,(j+ns)*(2*nv+1)+k+nv+1)+ &
            pf(i,5)*pf(i,6)*pf(i,7)*(strike*cnd(-d2)-f*cnd(-d1))-v0
        end if
    end do
  end do
end do

do i=1,n
  strike=lq(i,1)
  f=lq(i,3)
  lm=log(f/strike)
  sig=lq(i,8)+lq(i,9)*lm+lq(i,10)*lm**2+lq(i,11)*lm**3+lq(i,12)*lm**4+ &
      lq(i,13)*lm**5+lq(i,14)*abs(lm)
  d1=(log(f/strike)+(sig**2/2.)*lq(i,4))/(sig*sqrt(lq(i,4)))
  d2=d1-sig*sqrt(lq(i,4))
  if (lq(i,2)==0)then
    v0=lq(i,5)*lq(i,6)*lq(i,7)*(f*cnd(d1)-strike*cnd(d2))
  else
    v0=lq(i,5)*lq(i,6)*lq(i,7)*(strike*cnd(-d2)-f*cnd(-d1))
  end if
  do j=-ns,ns
    do k=-nv,nv
      f=lq(i,3)*(1+j*ds)
      lm=log(f/strike)
      sig=lq(i,8)+lq(i,9)*lm+lq(i,10)*lm**2+lq(i,11)*lm**3+lq(i,12)*lm**4+ &
          lq(i,13)*lm**5+lq(i,14)*abs(lm)
      sig=sig+lq(i,8)*k*dv
      d1=(log(f/strike)+(sig**2/2.)*lq(i,4))/(sig*sqrt(lq(i,4)))
      d2=d1-sig*sqrt(lq(i,4))
      if (lq(i,2)==0)then
        ml((j+ns)*(2*nv+1)+k+nv+1,i)=lq(i,5)*lq(i,6)*(f*cnd(d1)-strike*cnd(d2))-v0
      else
        ml((j+ns)*(2*nv+1)+k+nv+1,i)=lq(i,5)*lq(i,6)*(strike*cnd(-d2)-f*cnd(-d1))-v0
      end if
    end do
  end do
end do

f=lq(1,3)
do j=-ns,ns
  do k=-nv,nv
    ml((j+ns)*(2*nv+1)+k+nv+1,n+1)=lq(1,6)*f*j*ds
  end do
end do

call inverse(matmul(transpose(ml),ml),matinv,n+1)
mat_post=matmul(ml,matinv)
alloc=matmul(vit,mat_post)

end subroutine

double precision function cnd(x)
! standard normal distribution for computing BS formula
double precision ,parameter :: dpi=3.141592653589793238
double precision           :: x
double precision           :: l, k
double precision           :: a1,a2,a3,a4,a5
a1=0.31938153
a2=-0.356563782
a3=1.781477937
a4=-1.821255978
a5=1.330274429
l=abs(x)
k=1./(1.+0.2316419*l)
cnd=1.-1./Sqrt(2.*dpi)*Exp(-l**2./2.)*(a1*k+a2*k**2.+a3*k**3.+a4*k**4.+a5*k**5.)
If (x<=0) then
  cnd=1.-cnd
End If
end function

subroutine inverse(a,c,n)
!============================================================
! Inverse matrix
! Method: Based on Doolittle LU factorization for Ax=b
! Alex G. December 2009
!-----------------------------------------------------------
! input ...
! a(n,n) - array of coefficients for matrix A
! n      - dimension
! output ...
! c(n,n) - inverse matrix of A
! comments ...
! the original matrix a(n,n) will be destroyed
! during the calculation
!===========================================================
implicit none
integer n
double precision a(n,n), c(n,n)
double precision L(n,n), U(n,n), b(n), d(n), x(n)
double precision coeff
integer i, j, k

! step 0: initialization for matrices L and U and b
! Fortran 90/95 aloows such operations on matrices
L=0.0
U=0.0
b=0.0

! step 1: forward elimination
do k=1, n-1
  do i=k+1,n
    coeff=a(i,k)/a(k,k)
    L(i,k) = coeff
    do j=k+1,n
       a(i,j) = a(i,j)-coeff*a(k,j)
    end do
  end do
end do

! Step 2: prepare L and U matrices
! L matrix is a matrix of the elimination coefficient
! + the diagonal elements are 1.0
do i=1,n
  L(i,i) = 1.0
end do
! U matrix is the upper triangular part of A
do j=1,n
  do i=1,j
    U(i,j) = a(i,j)
  end do
end do

! Step 3: compute columns of the inverse matrix C
do k=1,n
  b(k)=1.0
  d(1) = b(1)
  ! Step 3a: Solve Ld=b using the forward substitution
  do i=2,n
    d(i)=b(i)
    do j=1,i-1
      d(i) = d(i) - L(i,j)*d(j)
    end do
  end do
  ! Step 3b: Solve Ux=d using the back substitution
  x(n)=d(n)/U(n,n)
  do i = n-1,1,-1
    x(i) = d(i)
    do j=n,i+1,-1
      x(i)=x(i)-U(i,j)*x(j)
    end do
    x(i) = x(i)/u(i,i)
  end do
  ! Step 3c: fill the solutions x(n) into column k of C
  do i=1,n
    c(i,k) = x(i)
  end do
  b(k)=0.0
end do
end subroutine inverse

Answer 1

这绝不是一个解决方案，但是 - 当您在 FORTRAN 中看到按 10e300 顺序排列的值（或无法重现的结果）时，通常意味着变量（或数组等）的值在) 未初始化。根据编译器设置等，FORTRAN 中声明的变量接收默认随机值（在gfortran 的情况下，顺序为10e300 或10e-300）。如果您要除以较小的（10e-300 或大约），这可能会给您NaN（或错误取决于编译器选项），因为10e-300 将被视为0 以达到工作精度。

我的建议：在您的FORTRAN 代码中添加一些打印语句并通过Python 调用它。确保在FORTRAN 子例程中声明的所有变量都已初始化。 （从涉及任何种类的划分开始）。