Interface legacy Fortran code from Python with f2py

I recently found some old Fortran code from the paper Algorithm AS197 A Fast Algorithm for the Exact Likelihood of Autoregressive-Moving Average Models (Melard 1984) [pdf] [code].

I’m not very fluent in Fortran77 and I don’t like the idea of rewriting a well tested piece of code introducing new bugs, so why not give a try to F2PY?

F2PY is a tool that provides an easy connection between Python and Fortran languages. F2PY is part of NumPy.

F2PY creates extension modules from (handwritten or F2PY generated) signature files or directly from Fortran sources.

The generated extension modules facilitate:

  • Calling Fortran 77/90/95, Fortran 90/95 module, and C functions from Python.
  • Accessing Fortran 77 COMMON blocks and Fortran 90/95 module data (including allocatable arrays) from Python.
  • Calling Python functions from Fortran or C (call-backs).
  • Automatically handling the difference in the data storage order of multi-dimensional Fortran and Numerical Python (i.e. C) arrays.

And this is all true!

First step, try to run f2py on the unmodified code:

$ f2py -c as197.f -m as197_

This command wraps and compile the Fortran code in as197.f and creates a module named as197_, without any hints the resulting module is importable and has a informative documentation:

flikam - Function signature:
  flikam(p,mp,q,mq,w,e,sumsq,fact,vw,vl,vk,toler,ifault,[n,mrp1,mr])
Required arguments:
  p : input rank-1 array('d') with bounds (1)
  mp : input int
  q : input rank-1 array('d') with bounds (1)
  mq : input int
  w : input rank-1 array('d') with bounds (n)
  e : input rank-1 array('d') with bounds (n)
  sumsq : input float
  fact : input float
  vw : input rank-1 array('d') with bounds (mrp1)
  vl : input rank-1 array('d') with bounds (mrp1)
  vk : input rank-1 array('d') with bounds (mr)
  toler : input float
  ifault : input int
Optional arguments:
  n := len(w) input int
  mrp1 := len(vw) input int
  mr := len(vk) input int

As you can see f2py has correctly wrapped most arrays making the length argument optional. Something strange is happening for p and q were the length is one. But looking at the Fortran code, that is the effect of a workaround for some VAX/VMS compiler.

C ORIGINAL CODE FROM STATLIB USED P(MP),Q(MQ) BELOW, BUT THIS
C FAILS ON COMPUTERS LIKE VAX/VMS WHEN MP OR MQ ARE ZERO.
      DIMENSION P(1),Q(1),W(N),E(N),VW(MRP1),VL(MRP1),VK(MR)

We can revert that part to the original one with the correct DIMENSION.

Another build and the output of:

python -c "import as197_; print as197_.flikam.__doc__"

is better than before:

flikam - Function signature:
  flikam(p,q,w,e,sumsq,fact,vw,vl,vk,toler,ifault,[mp,mq,n,mrp1,mr])
Required arguments:
  p : input rank-1 array('d') with bounds (mp)
  q : input rank-1 array('d') with bounds (mq)
  w : input rank-1 array('d') with bounds (n)
  e : input rank-1 array('d') with bounds (n)
  sumsq : input float
  fact : input float
  vw : input rank-1 array('d') with bounds (mrp1)
  vl : input rank-1 array('d') with bounds (mrp1)
  vk : input rank-1 array('d') with bounds (mr)
  toler : input float
  ifault : input int
Optional arguments:
  mp := len(p) input int
  mq := len(q) input int
  n := len(w) input int
  mrp1 := len(vw) input int
  mr := len(vk) input int

Now from the paper we know that some arguments are input, some output and some are workspace (so that no memory allocation is done inside the procedure).

The easy part is to use intent(in) or intent(out) for inputs and outputs, and intent(hide) for the array lengths. Then I expected to be forced to write a python wrapper to allocate all the workspace arrays but… f2py can do this for you!

Cf2py intent(in) P
Cf2py intent(hide) MP
Cf2py intent(in) Q
Cf2py intent(hide) MQ
Cf2py intent(in) W
Cf2py intent(out) E
Cf2py intent(hide) N
Cf2py intent(out) SUMSQ
Cf2py intent(out) FACT
Cf2py intent(hide,cache) VW
Cf2py intent(hide,cache) VL
Cf2py integer intent(hide) :: MRP1 = max(MP,MQ+1)+1
Cf2py intent(hide,cache) VK
Cf2py integer intent(hide) :: MR = max(MP,MQ+1)
Cf2py double precision optional :: TOLER=-1
Cf2py intent(out) IFAULT

We can specify intent(hide,cache) to the workspace arrays, set the computed dimension and f2py will generate all the bookkeeping code automatically.

We can check if the expanded signature is ok with:

f2py -h stdout as197.f

subroutine flikam(p,mp,q,mq,w,e,n,sumsq,fact,vw,vl,mrp1,vk,mr,toler,ifault) ! in as197.f
    double precision dimension(mp),intent(in) :: p
    integer optional,intent(hide),check(len(p)>=mp),depend(p) :: mp=len(p)
    double precision dimension(mq),intent(in) :: q
    integer optional,intent(hide),check(len(q)>=mq),depend(q) :: mq=len(q)
    double precision dimension(n),intent(in) :: w
    double precision dimension(n),intent(out),depend(n) :: e
    integer optional,intent(hide),check(len(w)>=n),depend(w) :: n=len(w)
    double precision intent(out) :: sumsq
    double precision intent(out) :: fact
    double precision dimension(mrp1),intent(hide,cache),depend(mrp1) :: vw
    double precision dimension(mrp1),intent(hide,cache),depend(mrp1) :: vl
    integer optional,intent(hide),depend(mq,mp) :: mrp1=max(mp,mq+1)+1
    double precision dimension(mr),intent(hide,cache),depend(mr) :: vk
    integer optional,intent(hide),depend(mq,mp) :: mr=max(mp,mq+1)
    double precision optional :: toler=-1
    integer intent(out) :: ifault
end subroutine flikam

Now with a minimal amount of coding, in fact we only declared a few f2py attributes, we have a working wrapper (the error management is still missing) with a nice signature:

flikam - Function signature:
  e,sumsq,fact,ifault = flikam(p,q,w,[toler])
Required arguments:
  p : input rank-1 array('d') with bounds (mp)
  q : input rank-1 array('d') with bounds (mq)
  w : input rank-1 array('d') with bounds (n)
Optional arguments:
  toler := -1 input float
Return objects:
  e : rank-1 array('d') with bounds (n)
  sumsq : float
  fact : float
  ifault : int

Here the full code listing.