Код Fortran дает неправильный результат при вызове из Python

Question

Чтобы улучшить скорость выполнения для финансовой задачи, я кодировал основные числовые части в Фортране, делая только доступ к файлам и т. Д. В Python. Я скомпилировал с помощью f2py, и вызвать подпрограмму fit (см. ниже в моем посте)

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

Я протестировал код на Фортране, в нем нет ошибок, однако после компиляции с помощью 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

Это ни в коем случае не решение, НО - когда вы видите значения в порядке 10e300 (или результаты, которые не воспроизводятся в этом отношении) в FORTRAN это обычно означает, что значения переменной (или массива и т. д.) не инициализируются. В зависимости от настроек компилятора и т. Д. Объявленная переменная в FORTRAN получает случайное значение по умолчанию (в случае gfortran в порядке либо 10e300 или же 10e-300). Если вы делите на меньшее (10e-300 или около того) это может дать вам NaN (или ошибка в зависимости от параметров компилятора) как 10e-300 будет рассматриваться как 0 с рабочей точностью.

Мой совет: поместите несколько печатных заявлений в свой FORTRAN код и позвонить через Python, Убедитесь, что все объявленные переменные в FORTRAN подпрограммы инициализируются. (Начните с тех, которые занимаются делением любого рода).