DOC: clarification about the interaction of numpy / f2py with Fortran procedures

[nbehrnd]

Issue with current documentation:

Prior to filing a PR about the presentation of f2py here, I would like to reach out for the maintainers of numpy for some discussion.

Limiting to the first example "the quick way", fib1.f currently presented is written in FORTRAN77 equally known as fixed style. Modern Fortran (everything since and including Fortran 90, with a last revision 2018, and standard 2023 in preparation) follows set of modernized rules, the statement implicit none might be one known equally by those outside Fortran. I rewrote the procedure as in snippet fib1a.f90 which f2py by

python -m numpy.f2py -c fib1a.f90 -m fib

accepts to wrap into an extension module. Python can use it to replicate replicate the computation a = np.zeros(8, 'd') to yield [0. 1. 1. 2. 3. 5. 8. 13.].

Observation 1: Line 6 defines double precision as suitable for the gfortran compiler which is applied on the elements of array a(n). With the advent of Fortran 2003, iso_fortran_env (ref) offers a compiler independent definition of precision which can be used instead. While fine for gfortran, and apparently fine for f2py, the subsequent use of the .so extension module by Python did not follow through with either version fib2a.f90 (which loads from iso_fortran_env only what is relevant to for "double precision" of a floating number), or fib2b.f90 (which does not impose such a constrain what other procedure definitions this Fortran intrinsic module might contain). For either one of the two, their use in the minimal Python example terminates with the report

[0. 0. 0. 0. 0. 0. 0. 0.]
corrupted size vs. prev_size
Aborted

Question: Does the use of iso_fortran_env in the Fortran procedure require a different syntax in the Python script to use the wrapper module? Does d in a = np.zeros(8, 'd') relate to double precision of reals (Fortran) / floating numbers (numpy)?

---

Observation 2: Especially Fortran subroutines are safer in use with an explicit definition if processed data shall

• only enter a procedure without being altered by this procedure (intent(in))
• shall return to the main program, possibly altered (intent(out)), or
• if the exchange between main program and routine shall be freely bidirectional, including the possibility that the procedure alters these values (intent(in out), or functional equivalent, intent(inout)). Hence I derived fib3a.f90 from fib1a.f90 which compiles well and yields the anticipated result [ 0. 1. 1. 2. 3. 5. 8. 13.]. However the use of these keywords and the double precision provided by Fortran's intrinsic iso_fortran_env module now yields the result [0. 0. 0. 0. 0. 0. 0. 0.]

Question: Here, I do not know how to interpret the result for the second case. Do you recommend a different input syntax on level of numpy?

---

The observations refer to an instance of Linux Debian 12/bookworm with Python (3.11.2) with numpy (1.24.2), gcc and gfortran (12.2.0-14) as provided by this Linux distributions. The Fortran modules were considered as "good enough" because their compilation to yield object files in a pattern of

gfortran ./fib3a.f90 -c -std=f2018 -Wall -Wextra -Wpedantic

did not yield an error. The two warnings e.g.

$ gfortran ./fib3a.f90 -c -std=f2018 -Wall -Wextra -Wpedantic
./fib3a.f90:16:15:

   10 |    do i = 1, n
      |              2 
......
   16 |       a(i) = a(i - 1) + a(i - 2)
      |               1
Warning: Array reference at (1) out of bounds (0 < 1) in loop beginning at (2) [-Wdo-subscript]
./fib3a.f90:16:26:

   10 |    do i = 1, n
      |              2            
......
   16 |       a(i) = a(i - 1) + a(i - 2)
      |                          1
Warning: Array reference at (1) out of bounds (-1 < 1) in loop beginning at (2) [-Wdo-subscript]

were noted, though not considered as relevant for the discussion. All relevant source code is in the zip archive below.

2023-04-21_f2py_discussion.zip

Idea or request for content:

No response

[HaoZeke]

Hi, sorry for getting to this belatedly. At the moment the intrinsic iso_fortran_env is not really supported, I recall adding some basic tests for compilation, will look into this.

[HaoZeke]

Ah, right so we currently 'support' the iso_fortran_env in the context of the type dictionary users can supply, via .f2py_f2cmap. To fix the compilation of fib2a.f90:

! file: fib2a.f90
subroutine fib(a, n)
   ! calculate first Fibonacci numbers
   use iso_fortran_env, only: dp => real64
   implicit none
   integer :: i, n
   real(kind=dp) :: a(n)

   do i = 1, n
   if (i .eq. 1) then
      a(i) = 0.0_dp
   elseif (i .eq. 2) then
      a(i) = 1.0_dp
   else
      a(i) = a(i - 1) + a(i - 2)
   end if
   end do
end subroutine fib
! end file fib2a.f90

We can write a type specification file .f2py_f2cmap given as:

dict(real=dict(dp='double'))

Note that this is kind of a hack, but it does work:

f2py -c fib2a.f90 -m fib2 --f2cmap .f2py_f2cmap
python -c "import numpy as np; import fib2; a = np.zeros(8, 'd'); fib2.fib(a); print(a)"
[ 0.  1.  1.  2.  3.  5.  8. 13.]

Other type specifications can be written out similarly
I guess this is then something to be clarified in the documentation.

Let me know if there are more questions about this, @nbehrnd .

[nbehrnd]

@HaoZeke Thank you for getting back in touch. I equally had to revisit the suggestion. Not knowing well enough the inner working of numpy your description of "kind of a hack, but it does work" is understood as you / as other maintainers of numpy would prefer to have some additional time to adjust how numpy internally works and the approach is a preliminary one.

This is why I returned to the other approach seen to explicitly define the numerical precision in Fortran independent to iso_fortran_env, i.e. in lines of

integer, parameter :: dp = selected_real_kind(15, 50)

as in fib1a.f90 perhaps easier for numpy to interact with. With the intent to substitute the snippets of old FORTRAN77 by more contemporary Fortran90+, I imagine already this could be more welcoming to the users of numpy. (It is influenced by Daniel Price' Fortran lectures with his recommendation "if there is a choice, prefer the modern Fortran".)

Both the original snippet of FORTRAN as well as the slightly modernized one however seem to interact with current numpy differently, than in the showcase of getting started/the quick way. The instructions in the info box successfully create a an array of eight zeros, to fill only the first six fields however fail. Because n no longer defaults to len(a) (as reported by print(fib1.fib.__doc__)), but to shape(a, 0) I speculate, could it be shape(a, 1) were a better (default) choice for a 1D array?

2023-05-22_f2py_discussion.zip