Procedural pointer in fortran

Question

Let us say I have the following abstract interface to a double precision function of single argument

module abstract

  abstract interface
     function dp_func (x)
       double precision, intent(in) :: x
       double precision :: dp_func
     end function dp_func
  end interface

end module abstract

In a different module I define two functions, a simple one g of the type dp_func and a more complicated one f

module fns

 contains
 double precision function f(a,b,x)
   double precision, intent(in)::a,b,x
   f=(a-b)*x 
 end function f

 double precision function g(x)
   double precision, intent(in)::x
   g=x**2 
 end function g
end module fns

Now a pointer to g can be created as follows

program main
use abstract,fns
procedure(dp_func), pointer :: p
double precision::x=1.0D0, myA=1.D2, myB=1.D1, y
p => g
y=p(x)
end program main

But how one can create a pointer to f(myA,myB,x), i.e., to f at fixed values of a and b, which can be regarded as a function of just 1 parameter, that is, of the dp_func type? At the end result I want to be able to write something like

p=>f(myA, myB, )
y=p(x)

Comments below suggest that function closure is not a part of fortran standard and that a wrapper function would be a possible solution to it. However, the wrapper must be initialized and this introduces some chances that end user may forget to call the initializer. How one can do it in a clean and transparent way?

EDIT After posting this question and googling with "closure and fortran", I found this example

which I present in picture form to emphasize the highlighting. This was presented in an online course. But I doubt such implicit parameter setting is a good programming practice. In fact, dangling variables like z in this example are perfect sources of errors!

Answer 1

You can use internal functions to wrap your functions, e.g.

program main
  use abstract
  use fns
  implicit none
  
  procedure(dp_func), pointer :: p
  double precision :: x, myA, myB, y
  
  x = 1.0D0
  myA = 1.D2
  myB = 1.D1
  
  p => g
  y=p(x)
  
  p => f2
  y = p(x) ! Calls f(1.D2, 1.D1, x)
  
  myA = 1.D3
  myB = 1.D2
  
  y = p(x) ! Calls f(1.D3, 1.D2, x)
contains
  double precision function f2(x)
    double precision, intent(in) :: x
    write(*,*) myA, myB
    f2 = f(myA,myB,x)
  end function
end program main

An internal function in a given scope can use variables from that scope, so they can act like closures.

The implicit use of myA and myB in the internal function f2 may well be a source of programming error, but, provided the scope of f2 is still in scope, this behaviour is identical to lambda functions in other languages, for example the equivalent python lambda:

f2 = lambda x: f(myA,myB,x)

As pointed out by @vladimirF, once the scope of f2 drops out of scope (e.g. if a pointer to f2 is stored and the procedure where f2 is declared returns) any pointers to f2 will become invalid. This can be seen in this code:

module bad
  use abstract
  use fns
  implicit none
contains

function bad_pointer() result(output)
  procedure(dp_func), pointer :: output
  
  double precision :: myA,myB
  
  myA = 1.D2
  myB = 1.D1
  
  output => f2
contains
  double precision function f2(x)
    double precision, intent(in) :: x
    write(*,*) myA, myB
    f2 = f(myA,myB,x)
  end function
end function

end module

program main
  use abstract
  use fns
  use bad
  implicit none
  
  procedure(dp_func), pointer :: p
  double precision :: y,x
  
  p => bad_pointer()
  x = 1.D0
  y = p(x)
end program

N.B. the above code may well run fine for this simple case, but it's relying on undefined behaviour so shouldn't be used.

Answer 2

You stated the following: "...However, the wrapper must be initialized and this introduces some chances that end user may forget to call the initializer. How one can do it in a clean and transparent way?..."

The following might be a solution. It still needs to be initialized but will throw errors if the user hasn't done so.

I defined a type closure which handles the function pointers.

! file closure.f90
module closure_m
  implicit none

  type closure
    private
    procedure(f1), pointer, nopass :: f1ptr => null()
    procedure(f3), pointer, nopass :: f3ptr => null()
    real :: a, b
  contains
    generic   :: init => closure_init_f1, closure_init_f3
      !! this way by calling obj%init one can call either of the two closure_init_fX procedures
    procedure :: exec => closure_exec
    procedure :: closure_init_f1, closure_init_f3
  end type

  abstract interface
    real function f1(x)
      real, intent(in) :: x
    end function

    real function f3(a, b, x)
      real, intent(in) :: a, b, x
    end function
  end interface

contains

  subroutine closure_init_f1(this, f)
    class(closure), intent(out) :: this
    procedure(f1)               :: f

    this%f1ptr => f
    this%f3ptr => null()
  end subroutine

  subroutine closure_init_f3(this, f, a, b)
    class(closure), intent(out) :: this
    procedure(f3)               :: f
    real,           intent(in)  :: a, b

    this%f1ptr => null()
    this%f3ptr => f
    this%a     =  a
    this%b     =  b
  end subroutine

  real function closure_exec(this, x) result(y)
    class(closure), intent(in) :: this
    real,           intent(in) :: x

    if      (associated(this%f1ptr)) then
      y = this%f1ptr(x)
    else if (associated(this%f3ptr)) then
      y = this%f3ptr(this%a, this%b, x)
    else
      error stop "Initialize the object (call init) before computing values (call exec)!"
    end if
  end function

end module

Concerning the lines class(closure), intent(out) :: this: This is the standard way of writing initializers for Fortran types. Note that it is class instead of type which makes this polymorphic as is needed for type-bound procedures.

I slightly adjusted your functions module (changed data types)

! file fns.f90
module fns_m
contains
  real function f(a, b, x)
    real, intent(in) :: a, b, x
    f = (a-b)*x
  end function

 real function g(x)
    real, intent(in) :: x
    g = x**2
 end function
end module

An example program

! file a.f90
program main
  use closure_m
  use fns_m

  implicit none

  type(closure) :: c1, c2

  call c1%init(g)
  print *, c1%exec(2.0)

  call c1%init(f, 1.0, 2.0)
  print *, c1%exec(2.0)

  call c2%init(f, 1.0, -2.0)
  print *, c2%exec(3.0)
end program

Example output

$ gfortran closure.f90 fns.f90 a.f90 && ./a.out
   4.00000000    
  -2.00000000    
   9.00000000