I am working on fractal generator code. The main code is written in python and iterations part is written in fortran. I use f2py to glue the two codes together. Here is the fortran function I use:
function iterate(z0, func, zmax, niter) result(n)
implicit none
complex(kind=8), intent(in) :: z0
real(kind=8), intent(in) :: zmax
integer, intent(in) :: niter
external func
complex(kind=8) :: func
integer :: n
complex(kind=8) :: z
n = 0
z = z0
do while ((n < niter) .and. (abs(z) < zmax))
z = func(z)
n = n + 1
end do
end function iterate
Here is the docstring for the generated wrapper code:
n = iterate(z0,func,zmax,niter,[func_extra_args])
Wrapper for ``iterate``.
Parameters
----------
z0 : input complex
func : call-back function
zmax : input float
niter : input int
Other Parameters
----------------
func_extra_args : input tuple, optional
Default: ()
Returns
-------
n : int
Notes
-----
Call-back functions::
def func(z): return z
Required arguments:
z : input complex
Return objects:
z : complex
I am getting Segmentation fault error when trying to
use iterate with any python callback function.
Here is a sample result that I get:
>>> from foo import iterate
>>> iterate(1.0j, lambda x: 4.0 + x**2, 4.0, 256)
Segmentation fault
I have looked through all available documentation on callbacks in f2py but haven't found any solution to this problem. Any help would be appreciated.
UPDATE
Here is a backtrace from gdb:
Program received signal SIGSEGV, Segmentation fault.
cb_func_in_iterate2__user__routines (return_value=0x7fffffffdbc0, z_cb_capi=0x3ff0000000000000)
at /tmp/tmpT8xG1q/src.linux-x86_64-2.7/juliamodule.c:470
470 /tmp/tmpT8xG1q/src.linux-x86_64-2.7/juliamodule.c: No such file or directory.
(gdb) backtrace
#0 cb_func_in_iterate2__user__routines (return_value=0x7fffffffdbc0, z_cb_capi=0x3ff0000000000000)
at /tmp/tmpT8xG1q/src.linux-x86_64-2.7/juliamodule.c:470
#1 0x00007ffff6b6482b in iterate2 (z0=(1,1), func=func@entry=0x7ffff6b60c20 <cb_func_in_iterate2__user__routines>, zmax=4, niter=256)
at julia.f90:38
#2 0x00007ffff6b64897 in f2pywrapiterate2 (iterate2f2pywrap=0, z0=(1,1), func=func@entry=0x7ffff6b60c20 <cb_func_in_iterate2__user__routines>,
zmax=4, niter=256) at /tmp/tmpT8xG1q/src.linux-x86_64-2.7/julia-f2pywrappers.f:25
#3 0x00007ffff6b61f5e in f2py_rout_julia_iterate2 (capi_self=<optimized out>, capi_args=<optimized out>, capi_keywds=<optimized out>,
f2py_func=0x7ffff6b64880 <f2pywrapiterate2>) at /tmp/tmpT8xG1q/src.linux-x86_64-2.7/juliamodule.c:811
#4 0x00000000004caaa1 in PyEval_EvalFrameEx ()
#5 0x00000000004c87a1 in PyEval_EvalCodeEx ()
#6 0x00000000005030ef in ?? ()
#7 0x00000000004f8c72 in PyRun_FileExFlags ()
#8 0x00000000004f7d77 in PyRun_SimpleFileExFlags ()
#9 0x00000000004982f2 in Py_Main ()
#10 0x00007ffff6f12b45 in __libc_start_main (main=0x497d80 <main>, argc=2, argv=0x7fffffffe2a8, init=<optimized out>, fini=<optimized out>,
rtld_fini=<optimized out>, stack_end=0x7fffffffe298) at libc-start.c:287
#11 0x0000000000497ca0 in _start ()
I've got a workaround, but not a complete fix, I'm afraid. The key point seems to be that f2py generates some strange looking, and possibly just incorrect, c code when func is a function. Changing func to be a subroutine seems to make it work.
To make your code work, change the following:
external func
complex(kind=8) :: func
...
do while ((n < niter) .and. (abs(z) < zmax))
z = func(z)
...
to
interface foo
subroutine func(f, z)
complex(kind=8), intent(out) :: f
complex(kind=8), intent(in) :: z
end subroutine func
end interface foo
...
do while ((n < niter) .and. (abs(z) < zmax))
call func(z, z)
...
This works fine for me, and gives the correct answer. If you're happy with that, you can stop reading. Below is how I worked out how to fix it.
As for why this happens, you need to run f2py like
f2py -m modname -h iterate.pyf iterate.f90
f2py -c --debug --build-dir build iterate.pyf iterate.f90
Then we can debug the code using gdb python, and then
run -c "import iterate; iterate.iterate(1.0j, lambda x: 4.0 + x*x, 4.0, 256)"
(although by the sounds of it, you probably got about this far anyway). Doing this, we find it segfaults on the following line in build/src.linux-x86_64.3.4/iteratemodule.c (your directory might be different):
complex_double z=(*z_cb_capi);
z_cb_capi being 0x0, hence the segfault. This appears in the Cython helper function that is actually called when func is called in the Fortran code. The function declaration appears just before the line above, and in the version where func is a function rather than a subroutine, it looks like the following mess:
static
#ifdef F2PY_CB_RETURNCOMPLEX
complex_double
#else
void
#endif
cb_func_in_iterate__user__routines (
#ifndef F2PY_CB_RETURNCOMPLEX
complex_double *return_value
#endif
#ifndef F2PY_CB_RETURNCOMPLEX
,
#endif
complex_double *z_cb_capi) {
which is essentially one of:
static void cb_func_in_iterate__user__routines (complex_double *return_value, complex_double *z_cb_capi)
static complex_double cb_func_in_iterate__user__routines (complex_double *z_cb_capi)
depending on whether or not F2PY_CB_RETURNCOMPLEX is defined. It appears that F2PY_CB_RETURNCOMPLEX never does get defined, because gdb reports the function as having the first definition, but it is being called as if it were the second. You can see this because return_value is set to the address of z in the Fortran code, and nothing (NULL) is passed to z_cb_capi.
So, that eventually leads us to the second way to fix this: pass -DF2PY_CB_RETURNCOMPLEX to f2py:
f2py -c -DF2PY_CB_RETURNCOMPLEX --debug --build-dir build iterate.pyf iterate.f90
This then compiles cb_func_in_iterate__user__routines in the second form, which means it does get called correctly.
Now, I don't think you should have to do this, so I rather suspect this is a bug in f2py.