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!

Fortran doesn't have _closures_, but if you can give an example of how you want to use this `f` then we can give specific examples of appropriate approaches. — francescalus, Apr 19 '21 at 11:44
@francescalus Thanks for a quick feedback. In a completely different module, which I do not want to modify, there is a subroutine that accepts only the `dp_func` pointers as a parameter. I want to provide such a pointer, but at the same time be flexible in the choice of `a` and `b`, that is, `a` and `b` are not known at the compile time. — yarchik, Apr 19 '21 at 11:51
You can likely extend the ideas of [this other question](https://stackoverflow.com/q/24127313/3157076) (and similar) with procedure pointers. — francescalus, Apr 19 '21 at 11:59
@francescalus Yes, that is the solution that I would typically use. I thought some new features have been added to new fortran standards that may enable closures. I am not comfortable with the idea of 1) writing a wrapper, 2) writing an initializer for a wrapper, 3) users forgetting to initialize before calling the wrapper. So in essence my question is how to indicate to the enduser that initializer must be called. — yarchik, Apr 19 '21 at 12:06
Fortran does not have closures. Period. You can create functor objects. https://stackoverflow.com/a/37744294/721644 — Vladimir F Героям слава, Apr 19 '21 at 20:13
@yarchik See my answer [below](https://stackoverflow.com/a/67168422/10774817). The object will throw an error if `exec` is called before initializing it. Does that help? — jack, Apr 21 '21 at 09:14

veryreverie · Answer 1 · 2021-04-21T11:03:59.703

4

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.

edited Apr 21 '21 at 11:03

answered Apr 19 '21 at 14:28

veryreverie

2,871
2
13
26

1

How does this approach differ from that highlighted (and rejected) in the question itself? – francescalus Apr 19 '21 at 15:06
@francescalus I started writing this answer before that edit to the question was made. – veryreverie Apr 19 '21 at 15:29
1

The most important point about a real closure is that it is a first class object that can be stored, passed around and so on. If you do pass around the pointer to the the internal function where the parent scope no longer exists, the context will be invalid, it is not stored in the pointer, unlike a real closure. – Vladimir F Героям слава Apr 19 '21 at 20:27
@VladimirF ah, thanks, I didn't know this. I should check this in my own code. – veryreverie Apr 19 '21 at 22:19
Can you please provide a small example where `f2` drops out of scope. – yarchik Apr 21 '21 at 10:29
@yarchik I have added an example which I believe gives an invalid pointer. – veryreverie Apr 21 '21 at 11:05
There is no way of having two procedure pointers using different parameters by using global variables `myA,myB`, is there? – jack Apr 21 '21 at 11:24
@jack I would be very wary of a solution involving global variables here. I think if a more robust solution were needed, your answer with a closure type would be a much better solution. – veryreverie Apr 21 '21 at 12:33

jack · Accepted Answer · 2021-04-21T11:45:19.907

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

Thank you jack. I noticed your answer and I like your approach. But I will wait a bit with accepting because I am curious to know if more people with come up with different solutions. As to your solution, can you briefly explain what is the meaning `class(closure), intent(out) :: this` and `generic :: init => closure_init_f1, closure_init_f3`. Yet another point, in the main program you define `type(closure) :: c`, can one define, say, `c1` and `c2`, which contain same functions, but store different context variables? — yarchik, Apr 21 '21 at 09:46
@yarchik I edited the answer and tried to include the questions from your comment. Feel free to comment again if it is still not self-explanatory. — jack, Apr 21 '21 at 11:16

Procedural pointer in fortran

2 Answers2

Linked