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

Have you tested being able to call a much simpler Fortran subroutine? — francescalus
– francescalus, Commented Sep 18, 2016 at 10:45
Polynomial evaluation by Horner method en.wikipedia.org/wiki/Horner%27s_method is less vulnerable to varying compiler optimizations and production of NaN/Inf — tim18
– tim18, Commented Sep 18, 2016 at 12:59
pf(i,8)+lm*((pf(i,9)+pf(i,10)*lm)+lm**2*(pf(i,11)+lm*(pf(i,12)+ & pf(i,13)*lm)) might be a better compromise — tim18
– tim18, Commented Sep 18, 2016 at 13:08
Fortran compilers are permitted to insert such grouping but don't optimize it. — tim18
– tim18, Commented Sep 18, 2016 at 13:11

ptev · Accepted Answer · 2016-09-18 10:53:32Z

2

This is by no means a solution, BUT - when you see values in the order 10e300 (or results that are not reproducible for that matter) in FORTRAN it usually means that the values of a variable (or array etc.) are not initialised. Depending on compiler settings etc. a declared variable in FORTRAN receives a default random value (in the case of gfortran in the order either 10e300 or 10e-300). If you are dividing by the smaller (10e-300 or thereabout) this may give you NaN (or an error depending on compiler options) as 10e-300 would be treated as 0 to working precision.

My advice: put some print statements in your FORTRAN code and call it with through Python. Make sure all declared variables in the FORTRAN subroutines are initialised. (Start from those which are involved in division of any kind).

answered Sep 18, 2016 at 10:53

ptev

1832 silver badges7 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

ZRH Over a year ago

ptev, thanks for your swift response, this is indeed a lead to look into! since my routine is time critical, i did not want to write code that loops over all array / vector elements. not being very experienced in fortran: is there a way i can set all elements of an array / vector to zero when i define it ?

ptev Over a year ago

Yup, fortran is probably the best language to do so in - it implicitly loops over arrays/vectors and does the operation to all members. Example A = 0.0 sets all the members of A to 0.0 (real zero). You can also do the value allocation when declaring the array, i.e. double precision :: vit(...) = 0.0. In this way also operations are done over entire arrays/vectors, i.e. sum(A) returns the sum of all elements of A. Good luck!

ZRH Over a year ago

Instead of outputting alloc, as in above code, i output vit and mat_post, which are absolutely reproducible in every run. however alloc, which is obtained via the last command line: alloc=matmul(vit,mat_post) is again not reproducible

ptev Over a year ago

Alright, so you have narrowed down the problem. Couple of things: When you declare the alloc I don't think you need to specify its size - it is an input and you should be able to leave it with a deferred shape :: alloc(:,:) . When using just Fortran you would also add implicit none after the subroutine line in its definition, but I don't know how this works with f2py. Finally, note that specifying intent (out) for alloc gets rid of its input values and keeps whatever is assigned in the subroutine. You can try removing that, not modifying alloc and check if input values are kept.

ZRH Over a year ago

issue solved now - i set all matrices to zero before running the computations & not seeing the abovementioned problems anymore. obtaining correct results reproducibly. i do not quite understand how, since intermediate results using the same arrays seemed reproducible, but anyway, issue solved ...

Collectives™ on Stack Overflow

Fortran code delivers wrong result when called from Python

1 Answer 1

5 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Related