Fortran code delivers wrong result when called from Python

Question

In order to improve speed of execution for a finance problem, I coded the core numerical parts in Fortran, doing only file access etc in Python. I compiled using f2py, and call the subroutine fit (see lower down in my post)

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

I have tested the Fortran code, there are no bugs in it, however, after compiling with f2py and calling it from Python, outputs are not reproducible. Sometimes, runs will come back with correct results (as compared to the identical algorithm programmed in Python), other times the call returns a vector with nan, sometimes results are on the order of 10^300. At present I have no idea how to identify the root cause, any help is highly appreciated !

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

Fortran code delivers wrong result when called from Python

Answers (1)

Related Questions