Reputation: 679
I'm trying to pass a python function to a f2py generated function. I can not figure out how to do this when the python function's input and output is an array. (I'm new to f2py!)
Here is my fortran code (num_Dg.f):
subroutine num_Dg(func,U,jacob,n)
double precision, intent(in) :: U(n)
double precision, intent(out) :: jacob(n,n)
double precision Up(n), Um(n)
double precision g_valp(n),g_valm(n)
double precision :: rel_eps = 2.e-5
double precision eps
integer i, j
external func
eps = minval(U)*rel_eps
do j=1,n
do i=1,n
if (i .eq. j) then
Up(i) = U(i)+(eps)
Um(i) = U(i)-(eps)
else
Up(i) = 0
Um(i) = 0
endif
enddo
call func(Up,g_valp,n)
call func(Um,g_valm,n)
do i=1,n
jacob(i,j) = 0.0
jacob(i,j) = (g_valp(i)-g_valm(i))/(2*eps)
enddo
enddo
end subroutine
I use f2py to compile: f2py -c num_Dg.f -m num_Dg
. It compiles without errors. But when import the function into python I get this:
>>> from num_Dg import num_dg
>>> print(num_dg.__doc__)
jacob = num_dg(func,u,[n,func_extra_args])
Wrapper for ``num_dg``.
Parameters
----------
func : call-back function
u : input rank-1 array('d') with bounds (n)
Other Parameters
----------------
func_extra_args : input tuple, optional
Default: ()
n : input int, optional
Default: len(u)
Returns
-------
jacob : rank-2 array('d') with bounds (n,n)
Notes
-----
Call-back functions::
def func(up,gvalp,[n]): return
Required arguments:
up : input rank-1 array('d') with bounds (n)
gvalp : input rank-1 array('d') with bounds (n)
Optional arguments:
n : input int, optional
Default: len(up)
There is a problem with the callback function func
. The input should be array U
, and the output should be gval
. Instead, both U
and gval
are inputs, with no output.
I'm assuming this is why the following python script returns garbage:
import num_Dg
def func(x):
return [x[0]**2,x[1],x[2]]
val = [0.2,0.1,0.1]
num_Dg.num_dg(func,val)
How should I change my fortran subroutine so that the callback function works properly? Any advice is appreciated! Thanks!
Upvotes: 1
Views: 359
Reputation: 679
I ended up figuring it out. I think the reason why it wasn't working was because I was missing the line cf2py intent(out), depend(n) :: g_valp
. The script can be compiled in the normal way with f2py: f2py -c num_Dg.f -m num_Dg
subroutine num_Dg(U,jacob,n)
cf2py intent(callback) func
external func
integer n
double precision :: U(n)
double precision :: jacob(n,n)
double precision Up(n), Um(n)
double precision g_valp(n),g_valm(n)
double precision :: rel_eps = 2.e-5
double precision eps
integer i, j
cf2py intent(in,copy), depend(n) :: U
cf2py intent(hide) :: n
cf2py intent(out), depend(n) :: jacob
cf2py intent(out), depend(n) :: g_valp
eps = minval(U)*rel_eps
do j=1,n
do i=1,n
if (i .eq. j) then
Up(i) = U(i)+(eps)
Um(i) = U(i)-(eps)
else
Up(i) = 0
Um(i) = 0
endif
enddo
call func(n,Up,g_valp)
call func(n,Um,g_valm)
do i=1,n
jacob(i,j) = (g_valp(i)-g_valm(i))/(2*eps)
enddo
enddo
end subroutine
Upvotes: 0