numpy - f2py, Python function that returns an array (vector-valued function)

Question

In the following Python I have five functions contained in the array returned by func which I have to integrate. The code calls an external Fortran module generated using f2py:

import numpy as np
from numpy import cos, sin , exp
from trapzdv import trapzdv
def func(x):
    return np.array([x**2, x**3, cos(x), sin(x), exp(x)])

if __name__ == '__main__':
    xs = np.linspace(0.,20.,100)
    ans =  trapzdv(func,xs,5)
    print 'from Fortran:', ans
    print 'exact:', np.array([20**3/3., 20**4/4., sin(20.), -cos(20.), exp(20.)])

The Fortran routine is:

      subroutine trapzdv(f,xs,nf,nxs,result)
          integer :: I
          double precision :: x1,x2
          integer, intent(in) :: nf, nxs
          double precision, dimension(nf) :: fx1,fx2
          double precision, intent(in), dimension(nxs) :: xs
          double precision, intent(out), dimension(nf) :: result
          external :: f 
          result = 0.0
          do I = 2,nxs
            x1 = xs(I-1)
            x2 = xs(I)
            fx1 = f(x1)
            fx2 = f(x2)
            result = result + (fx1+fx2)*(x2-x1)/2
          enddo
          return
      end 

The problem is that Fortran is only integrating the first function in func(x). See the print result:

from Fortran: [ 2666.80270721  2666.80270721  2666.80270721  2666.80270721  2666.80270721]
exact: [  2.66666667e+03   4.00000000e+04   9.12945251e-01  -4.08082062e-01 4.85165195e+08]

One way to workarond that is to modify func(x) to return the value of a given position in the array of functions:

def func(x,i):
    return np.array([x**2, x**3, cos(x), sin(x), exp(x)])[i-1]

And then change the Fortran routine to call the function with two parameters:

      subroutine trapzdv(f,xs,nf,nxs,result)
          integer :: I
          double precision :: x1,x2,fx1,fx2
          integer, intent(in) :: nf, nxs
          double precision, intent(in), dimension(nxs) :: xs
          double precision, intent(out), dimension(nf) :: result
          external :: f 
          result = 0.0
          do I = 2,nxs
            x1 = xs(I-1)
            x2 = xs(I)
            do J = 1,nf
                fx1 = f(x1,J)
                fx2 = f(x2,J)
                result(J) = result(J) + (fx1+fx2)*(x2-x1)/2
            enddo
          enddo
          return
      end 

Which works:

from Fortran: [  2.66680271e+03   4.00040812e+04   9.09838195e-01   5.89903440e-01 4.86814128e+08]
exact: [  2.66666667e+03   4.00000000e+04   9.12945251e-01  -4.08082062e-01 4.85165195e+08]

But here func is called 5 times more than necessary (in the real case func has above 300 functions, so it will be called 300 times more than necessary).

• Does anyone know a better solution to make Fortran recognizes all the array returned by func(x)? In other words, make Fortran build fx1 = f(x1) as an array with 5 elements corresponding to the functions in func(x).

---

OBS: I am compiling using f2py -c --compiler=mingw32 -m trapzdv trapzdv.f90

Answer 1

Unfortunately, you cannot return the array from the python function into Fortran. You would need a subroutine for that (meaning it is called with the call statement), and this is something that f2py does not let you do.

In Fortran 90 you can create functions that return arrays, but again this is not something that f2py can do, especially since your function is not a Fortran function.

Your only option is to use your looping workaround, or a redesign of how you want python and Fortran to interact.