Thread-safe initialization of shared module data (OpenMP)

Help
4

To the extent any engineer or scientist is pursuing programming with a compiler-based approach such as C++ or Fortran, especially Fortran given its legacy in Formula Translation, I always recommend to never lose sight of the underlying math and whether it can be described using constant expressions and to consequently pursue compile-time computing as much as possible.

Given the visionaries such as Stroustrup and Dos Reis constantly leading C++ forward, this is much easier and considerably more powerful presently with modern C++ than Fortran. But if one keeps seeking opportunities with Fortran and influences the language advancement with better facilities, more Fortran practitioners can benefit from the options.

So how is any of this relevant to the topic on hand? Well, it is to try to make the very question of shared data and whether it can thread-safe entirely moot, if possible. And an option to do so can be via named constants.
And note named constants really lend themselves well with MODULEs in Fortran.

In the case of the specific example in the original post, each entry, t_n, in the table is simply:

t_n = \frac{n!}{128^{n-1}}

which is a constant expression and thus can be encoded via an array named constant in Fortran. See the following code snippet that is spread over 4 lines mostly to make it easy on compilers, but also to make it easier for a future reader to follow the code. As to the former, the reason is the unfortunate situation where Fortran compilers do tend to struggle with constant expressions and throw ICEs, or worse (what can be worse than an ICE, right?), give wrong results:

      integer, parameter :: idx(*) = [( i, i = 1, n )]
      real(WP), parameter :: FAC = real(128, kind=WP)
      integer, parameter :: factorial(*) = [( product( idx(1:j) ), j = 1, n )]
      real(WP), parameter :: table(*) = [( real(factorial(i),kind=WP)/FAC**(i-1), i = 1, n )]

Of course I understand this may not always be the situation that program “data” can be characterized like so for thread-safe consumption. In such scenarios, one will have to pursue other techniques as described in this thread such as CRITICAL sections in legacy codes, or based on object-oriented (OO) design in modern codes. However, keeping an eye out for the named constant facility in Fortran is something I recommend.

Here’s a worked out example readers can try:

   integer, parameter :: WP = selected_real_kind( p=6 )
   integer, parameter :: n = 10
   integer :: i, j
   blk1: block
      integer, parameter :: idx(*) = [( i, i = 1, n )]
      real(WP), parameter :: FAC = real(128, kind=WP)
      integer, parameter :: factorial(*) = [( product( idx(1:j) ), j = 1, n )]
      real(WP), parameter :: table(*) = [( real(factorial(i),kind=WP)/FAC**(i-1), i = 1, n )]
      print *, "Block 1: array named constant"
      print *, "factorial up to ", n, ": ", factorial
      print *, "table: ", table
   end block blk1
   print *
   blk2: block
      real(WP) :: table(n)
      print *, "Block 2: using DO loop"
      call init_table( table )
      print *, "table: ", table
   end block blk2
contains
   subroutine init_table( t )
      real(WP), intent(inout) :: t(:)
      integer :: i
       t(1) = 1.0
       do i = 2, size(t)
         t(i) = real(i)*t(i-1)/real(128,kind=WP)
       end do
     end subroutine
end

Execution using gfortran:

C:\Temp>gfortran p.f90 -o p.exe

C:\Temp>p.exe
Block 1: array named constant
factorial up to 10 : 1 2 6 24 120 720 5040 40320 362880 3628800
table: 1.00000000 1.56250000E-02 3.66210938E-04 1.14440918E-05 4.47034836E-07 2.09547579E-08 1.14596332E-09 7.16227078E-11 5.03597164E-12 3.93435284E-13

Block 2: using DO loop
table: 1.00000000 1.56250000E-02 3.66210938E-04 1.14440918E-05 4.47034836E-07 2.09547579E-08 1.14596332E-09 7.16227078E-11 5.03597164E-12 3.93435284E-13

C:\Temp>

P.S.> Adapting the code for larger values of n and getting it to work with compilers’ of one’s choice is left up to the readers.

4 Likes