Here is a simple example that applies a pure (even elemental) function to each element of an array three different ways (elemental, do concurrent, openmp loop) and you get three different results:
program test_pure
implicit none
integer, parameter :: n = 16
integer :: i, A(n)
A = 1
A = h(A)
print *, A
A = 1
do concurrent (i = 1:n)
A(i) = h(A(i))
end do
print *, A
A = 1
!$OMP PARALLEL DO
do i = 1,n
A(i) = h(A(i))
end do
!$OMP END PARALLEL DO
print *, A
contains
elemental integer function h(x)
integer, intent(in) :: x
h = f()*x
end function
pure integer function f()
f = sum(A)
end function
end
This prints on my computer:
$ gfortran -fopenmp test_pure.f90
$ ./a.out
16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16
16 31 61 121 241 481 961 1921 3841 7681 15361 30721 61441 122881 245761 491521
241 481 61 121 92161 184321 16 31 23041 46081 961 1921 7681 7681 3841 368641
$ ./a.out
16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16
16 31 61 121 241 481 961 1921 3841 7681 15361 30721 61441 122881 245761 491521
7681 15361 16 31 61 121 481 961 1921 3841 7681 15361 53761 241 107521 215041
$ ./a.out
16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16
16 31 61 121 241 481 961 1921 3841 7681 15361 30721 61441 122881 245761 491521
245761 491521 61 121 961 1921 3841 7681 15361 30721 241 481 16 31 61441 122881
As you can see, the parallel loop even returns different results on every run!
See the thread What is a pure function? for a discussion of pure functions. A quick summary: functions can have two independent properties:
- deterministic: arguments fully determine the result (no reads from global variables unless they are compile time constants)
- side-effect-free: no side effects (no writes to global variables)
Fortran’s pure
is side-effect-free. Fortran’s 2023 simple
is both side-effect-free as well as deterministic.
It seems that currently “do concurrent” only requires side-effect-free, but not deterministic. As a consequence the above example returns non-deterministic results if run in parallel (I don’t have a compiler that can run do concurrent
in parallel, so I tested with openmp above). And even in serial you get different results depending on how you apply the elemental function. In particular the Fortran standard says about elemental functions:
In the array case, the values of the elements, if any, of the result are the same as would have been obtained if the scalar function had been applied separately, in array element order, to corresponding elements of each array actual argument.
This paragraph makes perfect intuitive sense to me. My understanding of it is that the elemental application above and the first loop should produce the same answer (on my computer a regular loop and do concurrent return the same answer). But that is not the case as you can see. If you look at the program, it seems to me it even cannot be the case. If the function is not deterministic, it seems in general you can’t guarantee the property that the result are the same as would have been obtained if the scalar function had been applied /.../ to corresponding elements of each array actual argument
and above is one such example, where the loop is not order independent, despite only using Fortran pure
functions. It seems one requires both side-effect-free as well as deterministic in order to ensure the loop is order independent.
What is the rationale to allow non-deterministic functions inside parallel loops and in elemental procedures?
Wouldn’t it make sense to only allow deterministic (and side-effect-free) functions?