How to pass parameters to a function in a thread-safe way?

aledinola · May 19, 2024, 9:14pm

Suppose I would like to integrate a function that depends on two parameters, alpha and beta. The function inputs should be x,alpha,beta, but the integrator (a numerical procedure that computes the integral) wants the function with one input only, x, as in here:

module integrals
implicit none
contains
function simpson(f, a, b) result(s)
! Taken from https://fortran-lang.org/en/learn/best_practices/callbacks/
  implicit none
  real(8), intent(in) :: a, b
  interface
  function f(x) result(func)
    real(8), intent(in) :: x
    real(8) :: func
  end function
  end interface
  real(8) :: s

  s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
end function simpson
end module integrals

One way of doing this is the following:

program main
use integrals, only: simpson
implicit none
real(8) :: alpha, beta
real(8) :: myint

alpha = 1d0
beta = 2d0

myint = simpson(myfunc, 0.0d0, 1.0d0)

contains

function myfunc(x) result(res)
    real(8), intent(in) :: x
    real(8) :: res
    res = alpha*x**2+beta
end function myfunc

end program main

So far so good, but the problem is that I would like to vary alpha and beta among several possible combinations, and to compute the integral of my function for each combo (alpha,beta). Something like this:

program main
use integrals, only: simpson
implicit none
integer, parameter :: n = 3
integer :: i1, i2
real(8) :: alpha, beta, alpha_vec(n), beta_vec(n), store(n, n)
real(8) :: myint

alpha_vec = [0.01d0, 2.0d0, 3.0d0]
beta_vec = [1.0d0, 0.0d0, 3.0d0]
!$omp parallel default(shared) private(i1,i2,myint,alpha,beta)
!$omp do collapse(2)
do i2=1, n
    do i1=1, n
        alpha = alpha_vec(i1)
        beta  = beta_vec(i2)
        myint     = simpson(myfunc, 0.0d0, 1.0d0)
        write(*,*) 'Integral = ', myint
        store(i1, i2) = myint
    enddo
enddo
!$omp enddo 
!$omp end parallel

write(*,*) 'Results:'
do i1=1, n
    write(*,*) store(i1, :)
enddo


contains

function myfunc(x) result(res)
    real(8), intent(in) :: x
    real(8) :: res
    res = alpha*x**2+beta
end function myfunc

end program main

Now, the code above works OK in serial, but as soon as I use openMP it gives garbage results. Reading online I found out that I should define alpha and beta as module variables, and declare them as threadprivate. Unfortunately I’m a bit lost, so I’d be very grateful if anyone has suggestions

Btw, I found very useful the “best practice” example on Callbacks — Fortran Programming Language but it does not talk about parallelization.

Beliavsky · May 19, 2024, 9:50pm

You could add a new derived type argument to simpson and f and put the parameters in it, for example

module common_mod
implicit none
type :: common_data
   real(8) :: c1, c2 ! extra parameters needed to evaluate f
end type
end module common_mod

module integrals
implicit none
contains
function simpson(f, a, b, dt) result(s)
! Taken from https://fortran-lang.org/en/learn/best_practices/callbacks/
  use common_mod, only: common_data
  real(8), intent(in) :: a, b
  type(common_data), intent(in) :: dt
  interface
  function f(x, dt) result(func)
    import common_data
    real(8), intent(in) :: x
    type(common_data), intent(in) :: dt
    real(8) :: func
  end function
  end interface
  real(8) :: s

  s = (b-a) / 6 * (f(a, dt) + 4*f((a+b)/2, dt) + f(b, dt))
end function simpson
end module integrals

aledinola · May 19, 2024, 9:59pm

Thanks! And then the derived type is declared as private?

Beliavsky · May 19, 2024, 10:14pm

No, because then the common_data derived type could not be used in modules other than common_mod.

FortranFan · May 19, 2024, 10:56pm

@aledinola ,

Refer to “Procedure Encapsulation Example” in Modern Fortran Explained by Metcalf et al. and also a section on thread-safe computation by @Arjen in “Modern Fortran in Practice”

Consider also employing a coarray Fortran program to independently unit test the code for thread-safety; that’s using MPI in a sense but with the far-simpler Fortran syntax.

Below is a silly example based on these concepts you can play around with. it is similar to the suggestion above re: the use of a derived type, only that it tries to make it a bit more generic and extensible,

Click to see example

module kinds_m
   use, intrinsic :: ieee_arithmetic, only : ieee_selected_real_kind
   integer, parameter :: WP = ieee_selected_real_kind ( p=12 )
end module
module user_func_m
   use kinds_m, only : WP
   private
   type, abstract, public :: user_func_t
      ! add members (components) as needed e.g., common numerical parameters
   contains
      procedure(IEval), pass(this), deferred :: Eval
   end type
   abstract interface
      elemental function IEval( this, x ) result( r )
         import :: WP, user_func_t
         class(user_func_t), intent(in) :: this
         real(WP), intent(in) :: x
         real(WP) :: r
      end function
   end interface
end module 
module integral_m
   use kinds_m, only : WP
   use user_func_m, only : user_func_t
   private
   type, abstract, public :: integral_t
      ! add members (components) as needed e.g., common numerical parameters
   contains
      procedure(IIntegrate), pass(this), deferred :: Integrate
   end type
   abstract interface
      elemental subroutine IIntegrate( this, uf, xlo, xhi, Integral )
         import :: WP, integral_t, user_func_t
         class(integral_t), intent(in) :: this
         class(user_func_t), intent(in) :: uf
         real(WP), intent(in) :: xlo
         real(WP), intent(in) :: xhi
         real(WP), intent(inout) :: Integral
      end subroutine
   end interface
end module 
module simpson_m
   use kinds_m, only : WP
   use user_func_m, only : user_func_t
   use integral_m, only : integral_t
   private
   type, extends(integral_t), public :: simpson_t
   contains
      procedure, pass(this) :: Integrate => Simpson_Integrate 
   end type
contains
   elemental subroutine Simpson_Integrate( this, uf, xlo, xhi, Integral )
      class(simpson_t), intent(in) :: this
      class(user_func_t), intent(in) :: uf
      real(WP), intent(in) :: xlo
      real(WP), intent(in) :: xhi
      real(WP), intent(inout) :: Integral
      Integral = (xhi - xlo) / 6 * ( uf%Eval(xlo) + 4*uf%Eval((xlo+xhi)/2) + uf%Eval(xhi) )
   end subroutine
end module 
module myfunc_m
   use kinds_m, only : WP
   use user_func_m, only : user_func_t
   private
   type, extends(user_func_t), public :: myfunc_t
      real(WP) :: alpha = 0.0_wp
      real(WP) :: beta = 0.0_wp
   contains
      procedure, pass(this) :: Eval => myfunc_Eval 
   end type
contains
   elemental function myfunc_Eval( this, x ) result( r )
      class(myfunc_t), intent(in) :: this
      real(WP), intent(in) :: x
      real(WP) :: r
      associate ( alpha => this%alpha, beta => this%beta ) 
         r = alpha*x**2 + beta
      end associate
   end function
end module
program p
   use kinds_m, only : WP
   use simpson_m, only : simpson_t
   use myfunc_m, only : myfunc_t
   type(simpson_t) :: simpson
   type(myfunc_t) :: myfunc[*]
   real(WP) :: Integral[*]
   integer :: i
   character(len=*), parameter :: fmtg = "(g0,t10,g0,t30,g0,t55,g0)"
   myfunc%alpha = 0.0_wp + 0.01_wp*this_image()
   myfunc%beta = 1.0_wp + 0.01_wp*this_image()
   call simpson%Integrate( myfunc, xlo=0.0_wp, xhi=1.0_wp, Integral=Integral )
   sync all
   if ( this_image() == 1 ) then
      ! Print results
      print fmtg, "Image", "Alpha", "Beta", "Integral" 
      do i = 1, num_images()
         print fmtg, i, myfunc[i]%alpha, myfunc[i]%beta, Integral[i] 
      end do
   end if 
end program

Program performance using the coarray unit test

C:\temp>ifx /free /standard-semantics /Qcoarray:shared p.f
Intel(R) Fortran Compiler for applications running on Intel(R) 64, Version 2024.1.0 Build 20240308
Copyright (C) 1985-2024 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.36.32537.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
Image    Alpha               Beta                     Integral
1        .1E-01              1.010000000000000        1.013333333333333
2        .2E-01              1.020000000000000        1.026666666666666
3        .3E-01              1.030000000000000        1.040000000000000
4        .4E-01              1.040000000000000        1.053333333333333
5        .5E-01              1.050000000000000        1.066666666666667
6        .6E-01              1.060000000000000        1.080000000000000
7        .7E-01              1.070000000000000        1.093333333333333
8        .8E-01              1.080000000000000        1.106666666666667
9        .9E-01              1.090000000000000        1.120000000000000
10       .1000000000000000   1.100000000000000        1.133333333333333
11       .1100000000000000   1.110000000000000        1.146666666666667
12       .1200000000000000   1.120000000000000        1.160000000000000

C:\temp>

egio · May 21, 2024, 12:19pm

An indirected call works also:

module integrals
implicit none
integer, parameter :: wp = kind(0.0d0) 

contains
    function simpson(f, a, b) result(s)
    ! Taken from https://fortran-lang.org/en/learn/best_practices/callbacks/
      implicit none
      real(wp), intent(in) :: a, b
      interface
      function f(x) result(func)
        import wp
        real(wp), intent(in) :: x
        real(wp) :: func
      end function
      end interface
      real(wp) :: s
    
      s = (b-a) / 6 * (f(a) + 4*f((a+b)/2) + f(b))
    end function simpson
end module integrals

module multi_integration
use integrals, only: simpson, wp
implicit none

contains

function multi(alpha, beta) result (s)
    
    real(wp) :: alpha, beta
    real(wp) :: s
    
    s = simpson(myfunc, 0.0d0, 1.0d0)
    
    contains
    
    function myfunc(x) result(res)
        real(wp), intent(in) :: x
        real(wp) :: res
        res = alpha*x**2+beta
    end function myfunc
end function

end module

program main
    use multi_integration
    implicit none
    integer, parameter :: n = 3
    integer :: i1, i2
    real(wp) :: alpha, beta, alpha_vec(n), beta_vec(n), store(n, n)
    real(wp) :: myint
    
    alpha_vec = [0.01d0, 2.0d0, 3.0d0]
    beta_vec = [1.0d0, 0.0d0, 3.0d0]
    !$omp parallel default(shared) private(i1,i2,myint,alpha,beta)
    !$omp do collapse(2)
    do i2=1, n
        do i1=1, n
            alpha = alpha_vec(i1)
            beta  = beta_vec(i2)
            myint     = multi(alpha, beta)
            write(*,*) 'Integral = ', i1, i2, myint
            store(i1, i2) = myint
        enddo
    enddo
    !$omp enddo 
    !$omp end parallel
    
    write(*,*) 'Results:'
    do i1=1, n
        write(*,*) store(i1, :)
    enddo
    
    end program main

aledinola · May 21, 2024, 11:25pm

Thanks! This works well and is very clear

Topic		Replies	Views
Defining an overloaded parameter? Help	4	184	January 30, 2025
Passing a subroutine as an argument to another sub Help	9	215	February 21, 2025
Best way to create a routine that can handle different precisions of input/output Help	11	875	April 5, 2023
Iso_c_binding: interface to a C function is not thread safe for openmp Help	2	661	May 19, 2021
Passing MIN to a function Help	4	372	June 15, 2023

How to pass parameters to a function in a thread-safe way?

Related topics