Implementation of a parametrized objective function without using module variables or internal subroutines

1

Below is a toy program that calls an optimization solver to optimize a simple parametrized objective function (see a similar example of root finding by @ivanpribec in a previous discussion if mine is not clear enough; see also the example in Note 12.18 on page 290 of WD 1539-1 J3/10-007r1, F2008 Working Document).

!----------------------------------------------------------------------!
! optim.f90: a simple program that calls an optimization solver to
! optimize a simple parameterized objective function.
!----------------------------------------------------------------------!


!----------------------------------------------------------------------!
! solver_mod: a module defining an optimization solver
module solver_mod

implicit none
private
public:: solver

! OBJ: an abstract interface for an objective function
abstract interface
    subroutine OBJ(x, y)
    real, intent(in):: x
    real, intent(out):: y
    end subroutine OBJ
end interface

contains

! solver: a doing-nothing solver for demonstration. It is provided by a 
! third library,  so  we CANNOT change its  signature to accept a hyper-
! parameter for `objective`.
subroutine solver(objective)
procedure(OBJ):: objective
real f
call objective(0.0, f)
end subroutine solver

end module solver_mod
!----------------------------------------------------------------------!


!----------------------------------------------------------------------!
program optimize

use solver_mod, only: solver
implicit none

! `hyper_parameter`: a hyperparameter for the objective function. In
! real applications, this parameter is not a constant, but a variable
! that cannot be predicted before the optimization starts.
real hyper_parameter
hyper_parameter = 42.0

call solver(objective)

contains

!--------------------------------------------------------------!
! objective: a simple objective function for demonstration.
! F2008 allows to pass internal procedures as actual arguments.
! See Note 12.18 on page 290 of WD 1539-1 J3/10-007r1 (F2008
! Working Document). We implement `objective` internally so that
! `hyper_parameter` is visible to it. We choose not to pass the
! parameter by a module variable, which is thread-unsafe.
subroutine objective(x, y)
real, intent(in):: x
real, intent(out):: y
y = (x+hyper_parameter)**2
end subroutine objective
!--------------------------------------------------------------!

end program optimize
!----------------------------------------------------------------------!

We assume that the code to call the solver is simply

call solver(objective)

In other words, solver accepts only one input, which is the subroutine that defines the objective function. In practice, it should also take inputs like the starting point and generate outputs such as the optimal point. We omit them for simplicity.

The signature of the objective function objective is simply

objective(x, f)

where x is the input, and f is the function value. However, the specific definition of objective depends on some hyperparameters that cannot be determined beforehand.

In our program, objective is implemented as an internal subroutine of the main program, so that the hyperparameter defined in the program is visible to objective. This enables us to pass the parameterized objective to solver without explicitly passing the hyperparameter. Fortran 2008 allows passing an internal subroutine to another subroutine as an actual argument. See Note 12.18 on page 290 of WD 1539-1 J3/10-007r1 (F2008 Working Document).

However, starting from gfrotran 13, such an internal subroutine will incur a linker warning about executable stack:

$ gfortran --version &&  gfortran optim.f90 -Wall

GNU Fortran (Ubuntu 13.3.0-6ubuntu2~24.04) 13.3.0
Copyright (C) 2023 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

/usr/bin/ld: warning: /tmp/cc4i8jW5.o: requires executable stack (because the .note.GNU-stack section is executable)

More seriously, it causes a segfault in MATLAB R2025a, which is probably a bug of MATLAB. MathWorks will take a long time to fix the bug if it is confirmed.

An alternative is to make the hyperparameter available to objective by a module variable, but this is not thread-safe recursion-safe.

Question: How to modify the implementation so that it is thread-safe without using internal subroutines without using module variables or internal subroutines? Of course, the implementation must be completely standard-conforming.

Any comments or criticism will be appreciated. Thank you very much.

The following is a very personal opinion. In terms of programming techniques, the example discussed here is almost the same as the one in Note 12.18 on page 290 of WD 1539-1 J3/10-007r1, F2008 Working Document. If an example in the Fortran standard leads to a safety issue (executable stack) in a major compiler (gfortran / Intel), this is something alarmingly wrong with the language or the compiler. The two cannot be simultaneously innocent in this case (sorry for complaining without providing a solution).

2 Likes
2

Not a direct response to your question, which requires careful reading, but having dealt with similar problems in the past, could execstack command resolve the issue? Here is a related discussion.

1 Like
3

You say that the value of hyperparameter is not known until the optimization starts. But does it change during the optimization?

If so, then hyperparameter is volatile and, regardless of whether the objective procedure is external or not, you must use some locking mechanism —e.g., pthread’s mutex or something like that.

Assuming no volatility is involved, the obvious thread-safe approach that does not depend upon internal procedures, is to use an abstract type, e.g.:

module solver_mod
    implicit none
    private

    type, abstract, public :: objective_t
    contains
        procedure(i_obj), deferred :: objective
    end type

    abstract interface
        subroutine i_obj(this, x, y)
            import
            class(objective_t), intent(in) :: this
            real, intent(in):: x
            real, intent(out):: y
        end subroutine
    end interface
    public :: i_obj

    public :: solver

contains
    subroutine solver(t)
        class(objective_t), intent(in) :: t
        real :: f
        call t%objective(0.0, f)
    end subroutine
end module solver_mod

module impl_mod
    use solver_mod

    implicit none
    private

    public :: solver_mod

    type, extends(objective_t), public :: impl
        real :: hyper_parameter
    contains
        procedure :: objective
    end type

contains
    subroutine objective(this, x, y)
        class(impl), intent(in) :: this
        real, intent(in):: x
        real, intent(out):: y
        print*,'optimizing...'
        y = (x+this%hyper_parameter)**2
    end subroutine objective
end module impl_mod

program optimize
    use impl_mod

    implicit none

    call solver(impl(12.))
end program optimize

The code above compiles properly with ifx, but you might have to adapt it to your compiler’s features.

3 Likes
4

Not sure if it will help you as a permanent solution, but if you’re using gfortran you could compile adding the flag -ftrampoline-impl=heap. That would avoid executable stack. In my case this was enough and it did not impact negatively in permormance (but maybe it could?)

4 Likes
5

No. It is a hyperparameter that defines the problem. It is not a parameter to be optimised.

6

Thank you @jwmwalrus . What is the signature of your solver? It must accept the objective function as the (only) input. Neither its signature nor that of the objective function can be modified by me, because they are defined by another library.

7

@jwmwalrus is suggesting the modern approach to your issue imho. Because your optimization problem is described (in general) by an objective function and a set of associated data, types, etc., best wrapping it with an abstract derived type as suggested above.

So your actual problem may actually look like:

type, public, extends(objective_t) :: my_problem
   real :: hyperparameter = 42.0
   contains
      procedure :: objective => my_objective
end type my_problem

...

        subroutine my_objective(this, x, y)
            class(my_problem), intent(in) :: this
            real, intent(in):: x
            real, intent(out):: y
            y = (x+this%hyper_parameter)**2
        end subroutine

and the solver interface will accept

subroutine solver(problem)
class(object_t), intent(inout) :: problem
call problem%objective(...) ! do stuff
end subroutine solver

This is always guaranteed 100% thread-safe.

3 Likes
8

Sorry, @FedericoPerini , but the interface of the solver is defined by another library and CANNOT be changed by me.

1 Like
9

Then maybe the problem lies in trying to achieve thread-related safety?

You mentioned that the actual code passes other arguments to solver. Are some of those passed to objective as well? Is any of those unique?

For the module-based approach, you’ll still need either some way to uniquely identify the thread in which objective is being executed, from arguments only; or to implement a queue (which implies locking).

11

No. The interface of objective is simply objective(x, f) and it cannot be changed.

To be more realistic, the interface of solver is solver(objective, x), where x inputs the starting point and outputs the solution. That’s all. It cannot be changed.

12

I agree that would be the best. I believe (but haven’t verified) that MacOS forbids having an executable stack, yet GCC compilers still work. I also read something about this in a recent C standard proposal: Functional Functions - A Comprehensive Proposal Overviewing Blocks, Nested Functions, and Lambdas for C,

NOTE: As of early 2025 in GCC 14, GCC selected a heap-based implementation that got rid of the executable stack. This requires some amount of dynamic allocation in cases where it cannot prove that the function is only passed down, not have its address taken in a meaningful way, or if it is not used immediately (as determined by the optimizer). It can be turned on with -ftrampoline-impl=heap.

@zaikunzhang, are you strictly limited to GCC 13?

Not all Fortran compilers use the executable stack approach. A brief discussion of the trampoline mechanism can be found here: Is creating nested subroutines/functions considered good practice in Fortran? - #25 by szakharin

Another discussion of issues related to internal procedures is from @sblionel in Doctor Fortran in “Think, Thank, Thunk”,

Historical note: in the past, these thunks were built on the stack, which was very convenient. However the notion of executable code on the stack was also convenient for virus writers so operating systems evolved to disallow this. Other solutions were found, a topic for someone else’s blog.

2 Likes
13

No, I am not limited to GCC 13. On the contrary, I aim to support ALL compilers. In other words, an issue with ANY compiler is not tolerable, let alone a mainstream one like GCC13.

14

Why using module variables is not thread safe? One can delcare them as threadprivate

I have this example on my github:

I have to minimize the function rhs using a minimization routine known as golden section search

subroutine sub_vfi()
use mod_param
use mod_numerical, only: golden_method
use omp_lib
implicit none
! Omitted code
call golden_method(rhs, small, c_ub-small, optim, f_max,1.0d-5)
  ! Omitted code
contains
    
    pure function rhs(c) 
        !Evaluates the RHS of the Bellman equation
        !This function could be put instead as a module function
        !in the module mod_baselib.f90
        implicit none
        ! Input arguments:
        real(8), intent(in) :: c
        ! Function result:
        real(8) :: rhs
        !Local variables:
        real(8) :: ap_next, v_next
        !Next-period assets
        ap_next = (1+r)*a_today+z_today-c
        v_next  = linint(asset_grid,EVz,ap_next)
        rhs = utilfun(c) + beta*v_next
   end function
end subroutine sub_vfi

The function rhs is defined as an internal function. Besides the variable to be chosen (the argmin) it has three variables z_today,a_today,EVz that are defined as globals in the module mod_param with the threadprivate attribute:

module mod_param
	implicit none
real(8), parameter :: large_negative = -1d100
!=================================================
  ! OMP DIRECTIVES
  !=================================================
  !stuff to be available to (sub)procedures and which may need to be varied across OpenMP threads
  !$omp threadprivate(z_today,a_today,EVz)
end module mod_param

I think Arjen Markus has a chapter dedicated to this problem in his nice book “Modern Fortran in Practice”

1 Like
15

I have recently ecountered a similar problem.
I discussed the issue with Gemini and it seemed to me that its final answer might be worth exploring. However, I confess I haven’t had the time to dig deeper.

3 Likes
16

I don’t think anyone has mentioned this approach yet in this discussion. You can define a shared array whose length is the number of OpenMP threads, and you can put things that are local to individual threads in that array. The function omp_get_num_threads() returns the number of threads. If it is a single integer, for example, then just make an integer array. If there are several variables, then make a derived type that collects all of them together. It is in shared memory, so it does not need to be passed as in argument. Use the omp_get_thread_num() value to index the correct element of the array (remember that the thread range is 0:N-1).

17

I tried the final answer by Gemini without checking it. The code does not compile.

A few other LLMs have also been tried.

  1. The code by Deepseek and ChatGPT 4 does not compile.
  2. Grok-4 first provides an implementation with a module variable, and then generates code that does not compile when requested to avoid module variables.
  3. Very interestingly, ChatGPT o3-pro claims that “there is no way in strictly standard-conforming Fortran to let objective see a run-time hyper-parameter and at the same time make that parameter private to every POSIX / OpenMP / TBB thread”.

See my prompt and ChatGPT o3-pro’s answer if you are interested. If it is right, I am sorry, this is a non-negligible limitation & flaw of Fortran as a scientific-computing language. Can anyone prove it wrong?

1 Like
18

In a previous reply, I asked about a unique way to identify the thread, because I kept thinking in terms of retrieving data from a Go-like goroutine-safe map —but I guess an array of channels works as well :slightly_smiling_face: .

19

@zaikunzhang I am trying to understand your question but I am not sure I understand it. To provide a callback in a thread safe manner, the standard way is to pass a context struct, and the external solver must pass it through the call chain into your callback.

You can also use an internal procedure. I tried GFortran 13 and 15 on macOS/arm, and I don’t get any warnings regarding executable stack, but not sure how they do it.

Besides the two approaches above, how else would you like to do it?

1 Like
20

Hi @certik . Maybe I have not emphasised the following point enough:

The solver is provided by another library that cannot be modified by me. The interface of the solver is fixed as solver(objective), and that of the objective is fixed as objective(x, f).

The question is essentially: what if I have an objective whose signature is bla(x, f, hyper_parameter), where hyper_parameter is unknown before the program starts? In other words, what if the implementation yof the objective depends on hyperparameters that are unknown beforehand? How to optimise such an objective using solver?

See my code above for details. See a similar example of root finding by @ivanpribec in a previous discussion if mine is not clear enough; see also the example in Note 12.18 on page 290 of WD 1539-1 J3/10-007r1, F2008 Working Document. As an example mentioned in the standard, it is not a wired situation but a common use case.

My approach is to use an internal procedure as you mentioned. But this leads to executable stack, which is a security issue. It is not allowed on some systems, and it causes a segfault of MATLAB R2025a.

Another approach is to pass the hyperparameter using a module variable, but it is not thread-safe or recursion-safe?

Is it possible to have a standard-confirming implementation that is thread-safe without internal procedures?

21

There is not a single occurence of the word “thread” in the Fortran standard (J3/24-007). So thread-safety is something beyond the Fortran standard.

Depending on the threading model (or library) in use, users can force make their program thread-safe, if they are aware the underyling library isn’t. This will most likely serialize the execution of their program, if the original purpose of parallelism was data-driven.

With OpenMP:

!$omp critical
hyper_parameter = 42.0
call solver(objective)
!$omp end critical

With pthreads using fortran-unix interfaces:

  type(c_pthread_mutex_t), save :: mutex
  integer :: stat

  ! only allow one thread to change hyper_parameter and objective callback at a time
  stat = c_pthread_mutex_lock(mutex)

  hyper_parameter = 42.0
  call objective(x,f)

  stat = c_pthread_mutex_unlock(mutex)

With TBB the pattern is similar (https://www.intel.com/content/www/us/en/docs/advisor/user-guide/2025-0/threading-building-blocks-simple-mutex-example.html).