Allow subroutine to accept any kinds of real variable without maintaining multiple subroutines

Help
1

I want to call Fortran subroutines in my python program through f2py. I know that f2py will implicitly make a copy of input array if the type of input array is inconsistent to the one declared in Fortran, and it will affect the performance.
I want to avoid the implicit copying, so I create two almost identical subroutines, one for single precision and another for double precision.
Here is an example of my code:

! sub.f90
! compiled by: f2py -c sub.f90 -m fortmod

subroutine do_something_r32(arr1, arr2, arrout)
    implicit none
    real*4, intent(in) :: arr1, arr2
    real*4, intent(out) :: arrout
    ...
    return
end subroutine do_something_r32

subroutine do_something_r64(arr1, arr2, arrout)
    implicit none
    real*8, intent(in) :: arr1, arr2
    real*8, intent(out) :: arrout
    ...
    return
end subroutine do_something_r64

# test.py
import numpy as np
from fortmod import do_something_r32, do_something_r64

def wrapper(arr1, arr2):
    dtype = arr1.dtype
    if dtype == np.float32:
        return do_something_r32(arr1, arr2)
    elif dtype == np.float64:
        return do_something_r64(arr1, arr2)

It works well.
The problem is, I have to maintain the two almost identical subroutines do_something_r32 and do_something_r64. They are almost identical, only the types of variables are different. So if I update the content of do_something_r32, I have to remember to update do_something_r64. It’s very inconvenience.
So I want to know that is there any method to deal with this problem? I hope that I only need to create and maintain one subroutine, but it can accept the different types of variable.
Thanks!

p.s.
I think that the most easy method is to transform the dtype of arr1 and arr2 in wrapper to the correct type before calling Fortran subroutine. But I found that it still cause some overhead so I want to deal with this problem at the Fortran end.

2

Use an INCLUDE statement in place of your …

Real*4 and Real*8 are not standard Fortran, they are a commonly found extension.

3

In your program, will do_something be invoked by both 32-bit and 64-bit floats, or only by one of the two? In general, are you using both 32-bit and 64-bit floats in the same program? In my programs I use a single type of float, defined in a module, so my code looks like this:

module kind_mod
use, intrinsic :: iso_fortran_env
implicit none
private
public :: wp
integer, parameter :: wp = real64 ! or real32 or real128
end module kind_mod

module sd_mod
use kind_mod, only: wp
implicit none
private
public :: sd_fixed_mean
contains
function sd_fixed_mean(xx,xmean) result(xsd)
! compute standard deviation around a fixed mean xmean (zero by default)
real(kind=wp), intent(in)           :: xx(:)
real(kind=wp), intent(in), optional :: xmean
real(kind=wp)                       :: xsd
real(kind=wp)                       :: xm
integer                             :: n
n = size(xx)
if (n < 1) return
if (present(xmean)) then
   xm = xmean
else
   xm = 0.0_wp
end if
xsd = sqrt(sum((xx-xm)**2)/n)
end function sd_fixed_mean
end module sd_mod
1 Like
4

Here is a proposal that would fix this particular issue:

But that is not implemented yet.

Until then, the best options are to use macros and/or includes. See this discussion in stdlib: How to implement same procedures for different numeric kinds · Issue #35 · fortran-lang/stdlib · GitHub where we decided to use the fypp macro processor.

5

Thanks for the helps, I use include statement suggested by @themos and it works well.
The first link provided by @certik is interesting! Hope that it can be included in fortran standard in the future. Or maybe include something similar to template in C++.

BTW, this is first time I learned that real*4 and real*8 are not standard Fortran… I learned Fortran in college several years ago and real*4 is what professor taught us. I have to say that an active forum is really helpful, avoiding me keep bad programming habit.

6

Before the advent of Fortran 90 (some thrity years ago ;)) the only way to define the precision was via the REAL and DOUBLE PRECISION keywords. The REAL4 and REAL8 declarations were an extension to these keywords and this extension was widely used. But as Themos stated, it has always been an extension.
With Fortran 90 the language got the KIND declaration, which is much more flexible. Note that the actual values that represent the various kinds are compiler-dependent. So to be robust, avoid hardcoded numbers like 4 and 8 (you often see that, sigh). Use the KIND() intrinsic or the SELECTED_*_KIND() intrinsic or the ISO_FORTRAN_ENV module instead.

7

What I know is that complex*8 and complex(8) are different.:laughing:

8

See this blog post by Steve Lionel: Doctor Fortran in "It Takes All KINDs" - Doctor Fortran

1 Like