Hi @ivanpribec , you have provided a perfect example supporting @rwmsu 's viewpoint quoted below, with which I totally agree.
1 Like
I don’t know enough compiler stuff to interpret output from compiler explorer, but @ivanpribec does, so I will bother him now
Thanks for the suggestion. I added a forall variant.
gfortran:
Variant Avg Time(s) Speedup (x) GFLOPS
Classic Nested Loops 0.8813E-06 1.000 22.7
Column-wise (1-loop) 0.8953E-06 0.984 22.3
Do Concurrent 0.8953E-06 0.984 22.3
...
Forall 0.8953E-06 0.984 22.3
flang:
Classic Nested Loops 0.4938E-06 1.000 40.5
Column-wise (1-loop) 0.5359E-06 0.921 37.3
Do Concurrent 0.4852E-06 1.018 41.2
...
Forall 0.1696E-05 0.291 11.8
One thing which I believe my example shows is that even a simple operation can be expressed in many different ways. It is a challenging task for the compiler to “see through” the programming language abstractions and deliver the optimally performing code. Not to mention that different optimizations may be needed depending on the sizes (m, n) (which may only be determinable at runtime) and for different hardware.
Procedures have always offered a convenient way of capturing the important building blocks, which can be optimized separately (or even implemented by other means than Fortran).
I’ve done a few experiments (and so has @hkvzjal) - if your elemental procedure is not in the same translation unit as the site where it is called, the compilers I’ve checked (gfortran, flang) fallback to calling the default scalar variant in a loop.
Here’s an example:
! stencil_helpers.f90
module stencil_helpers
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
elemental function stencil_point(c, w, e, s, n, cx, cy) result(val)
real(dp), intent(in) :: c, w, e, s, n, cx, cy
real(dp) :: val
val = c + cx*(e - 2.0_dp*c + w) + cy*(n - 2.0_dp*c + s)
end function stencil_point
end module
When compiled with gfortran -O3 -mcpu=native -S stencil_helpers.f90, the assembly contains a scalar function (we can infer it is scalar, because all registers start with d for double precision scalar):
___stencil_helpers_MOD_stencil_point:
LFB0:
ldr d31, [x0]
fmov d28, 2.0e+0
ldr d29, [x2]
ldr d30, [x4]
ldr d26, [x1]
fmsub d29, d31, d28, d29
ldr d27, [x3]
fmsub d30, d31, d28, d30
ldr d28, [x5]
ldr d0, [x6]
fadd d29, d29, d26
fadd d30, d30, d27
fmadd d31, d28, d29, d31
fmadd d0, d0, d30, d31
ret
In a second file we can now call the elemental function:
! stencil.f90
module stencil
use stencil_helpers, only: dp, stencil_point
implicit none
contains
subroutine stencil_elemental(uold, unew, nx, ny, coeffx, coeffy)
integer, intent(in) :: nx, ny
real(dp), intent(in) :: uold(0:nx+1,0:ny+1)
real(dp), intent(out) :: unew(0:nx+1,0:ny+1)
real(dp), intent(in) :: coeffx, coeffy
unew(1:nx,1:ny) = stencil_point( &
c=uold(1:nx,1:ny), &
w=uold(0:nx-1,1:ny), &
e=uold(2:nx+1,1:ny), &
s=uold(1:nx,0:ny-1), &
n=uold(1:nx,2:ny+1), &
cx=coeffx, cy=coeffy)
end subroutine stencil_elemental
end module
The assembly now contains a loop applying the scalar function, element after element:
___stencil_MOD_stencil_elemental:
; ... omitted ...
L3:
mov x0, x27
sub x1, x27, #8
add x27, x27, 8
add x4, x0, x19
add x3, x0, x20
mov x2, x27
add x6, x29, 152
add x5, x29, 144
bl ___stencil_helpers_MOD_stencil_point
str d0, [x26], 8
cmp x21, x27
bne L3
If we pack the two modules together and compile it (Compiler Explorer):
! stencil_fused.f90
module stencil
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
elemental function stencil_point(c, w, e, s, n, cx, cy) result(val)
real(dp), intent(in) :: c, w, e, s, n, cx, cy
real(dp) :: val
val = c + cx*(e - 2.0_dp*c + w) + cy*(n - 2.0_dp*c + s)
end function stencil_point
subroutine stencil_elemental(uold, unew, nx, ny, coeffx, coeffy)
use stencil_helpers, only: dp, stencil_point
integer, intent(in) :: nx, ny
real(dp), intent(in) :: uold(0:nx+1,0:ny+1)
real(dp), intent(out) :: unew(0:nx+1,0:ny+1)
real(dp), intent(in) :: coeffx, coeffy
unew(1:nx,1:ny) = stencil_point( &
c=uold(1:nx,1:ny), &
w=uold(0:nx-1,1:ny), &
e=uold(2:nx+1,1:ny), &
s=uold(1:nx,0:ny-1), &
n=uold(1:nx,2:ny+1), &
cx=coeffx, cy=coeffy)
end subroutine stencil_elemental
end module
the elemental function gets inlined and the compiler generates neat looking SIMD instructions (the v stands for vector registers, and postfix .2d means there are two doubles in the vector registers). Here it looks like it applied loop unrolling to a depth of 4:
.L6:
add x7, x7, 2
add x8, x8, 32
ldp q25, q24, [x6], 32
ldp q5, q2, [x2], 32
ldp q1, q23, [x5], 32
ldp q21, q26, [x4], 32
mov v4.16b, v25.16b
mov v0.16b, v24.16b
mov v3.16b, v5.16b
mov v20.16b, v2.16b
fmls v4.2d, v5.2d, v29.2d
fmls v0.2d, v2.2d, v29.2d
fmls v23.2d, v2.2d, v29.2d
fmls v1.2d, v5.2d, v29.2d
fadd v4.2d, v4.2d, v27.2d
fadd v0.2d, v0.2d, v25.2d
mov v27.16b, v24.16b
fadd v26.2d, v23.2d, v26.2d
fadd v21.2d, v1.2d, v21.2d
fmla v20.2d, v0.2d, v31.d[0]
fmla v3.2d, v4.2d, v31.d[0]
fmla v20.2d, v26.2d, v30.d[0]
fmla v3.2d, v21.2d, v30.d[0]
stp q3, q20, [x3], 32
cmp x9, x7
bne .L6
TL;DR, just because a function is marked as elemental, doesn’t mean it uses parallel evaluation.
9 Likes
very nice and informative talk @rouson (not surprised by it, of course). For those who missed it, here is a recording of the talk: https://www.youtube.com/watch?v=DKV2Whf4MKg
2 Likes
@ivanpribec thanks for all the great information on what’s being done with elemental. I’m especially glad for the tip about getting vectorization from putting the procedure and its invocation in one file. I wonder whether link-time optimization can help when the two are in different files. Also, I wonder if it would be reasonable to request a flag that multithreads or offloads elemental procedure invocations much like several compilers offer flags for multithreading or offloading do concurrent.
I suspect that a lot of what we face are chicken-and-egg problems. My sense is that compilers devote a lot of energy to optimizing loops because humans write lots of loops, but I bet more people would use elemental if compilers did better things with elemental. Once a programming practice becomes entrenched (initially for good reason because it’s the best option at the time), it can be very challenging to get compiler support for newer, less-entrenched practices.
3 Likes
but you can do
do concurrent(i=1:5,j=2:18,k=5:7,k<j)
elemental_fun(i,j,k)
end do
which might not be one line but it is not:
do i = 1,5
do j = 2, 18
do k = 5, 7
if k < j then
elemental_fun(i,j,k)
end if
end do
end do
end do
2 Likes
Regarding -flto/-ipo, yes it does help but one should be wary of two things: the actual generated ASM code might not be exactly the same between embedding the elemental function in the file (say with an #include or duplicating it) or using link time optimization to inline it from a different file unit. This might not be a deal breaker though. The actual deal breaker is that, using -flto/-ipo with real applications becomes practically almost impossible because build time becomes restrictive and the amount of build error one encounters can be overwhelming.
I personally don’t think that a flag to offload elemental functions would be of any good for practitioners, it would just bloat the already bloated ocean of obscure options. I personally believe that the most sane path forward for the regular practitioner would be that fortran compilers treated elementals as “inline”. Why don’t I think a flag would be a good path forward?: if an elemental procedure is bound to apply a scalar operation on all elements across all ranks of arrays, then it is the parent scope that shall pilot where the scalar operation is to be applied, meaning, if the operation is on the CPU the operation shall be done there, and if the parent scope is offloaded, the scalar function shall follow on the GPU. And this is basically because as of today, moving data between the GPU and the CPU is still very expensive for most of the available hardware, so if one is not careful about it, it can kill all performance. Fine grained control of data-movement is still needed.
2 Likes
This is a great discussion. I’m learning from several great insights about compilers.
FWIW, I’ll add a couple of philosophical points:
In general, I like features like elemental, do concurrent, and pure for things that don’t necessarily relate to with performance. I like what these features communicate to the reader. For example, pure corresponds to about 1.5 pages of constraints in the standard. That’s a lot of information about what the code can and can’t do. I find that information invaluable when I have to refactor code – especially if it’s code with which I’m not familiar. Things like code movement become a lot easier when I know, for example, that a function is pure and therefore cannot modify its arguments. Best of all is when I know that most or all functions in a library are pure because then I can move really fast when modifying code. To wit, all intrinsic functions are pure so we’re all used to using pure functions regularly and it’s nice when user-defined functions behavior similarly to intrinsic functions. To wit, the final sentence of Note 4 in Fortran 2023 Clause 15.7 “Pure procedures” states, “It is expected that most library procedures will conform to the constraints required of pure procedures, and so can be declared pure and referenced in DO CONCURRENT constructs, FORALL statements and constructs, and within user-defined pure procedures.”
I also like the that the use of these features can rule out certain mistakes. If I invoke an elemental procedure on an array, I know that I can’t make mistakes related to array bounds.
For the above reasons, I generally advocate writing most code first for clarity and robustness. Then profile the code and address bottlenecks, which might involve removing some of the aforementioned features in the few places that are truly performance-critical.
5 Likes
Talking of performance, I have made some progress with the GPU experiments I mentioned in this thread. But I have also encountered some disappointments 
. I will create a new thread for this.
1 Like
Following on from your talk Jane and I reran one of the examples from our fourth edition with the current crop of compilers. The example compares the times of evaluating a = b + c for a, b and c being one d comformant arrays. The example times the following constructs - simple do loop, whole array syntax, do concurrent and openmp parallel do. Here are the results for the Cray, gfortran, Intel ifx, Intel ifort, Nag and nvidia compilers.
| Fortran | Compiler | Cray | gfortran | Intel | Intel | Nag | nvidia |
| construct | ifort | ifx | |||||
| Version | 16.0.1 | 15.2.0 | 2021.13.0 | 2025.2.0 | 7.2.7225 | 25.11-0 | |
| Whole array | 2.0243 | 2.2956 | 1.9194 | 1.9312 | 4.2779 | 1.9447 | |
| Do loop | 2.5652 | 2.2878 | 1.8784 | 1.9203 | 2.0193 | 1.9312 | |
| do concurrent | 3.2638 | 2.2410 | 0.4900 | 0.5010 | 1.8422 | 1.9738 | |
| openmp | 0.9984 | 0.5010 | 0.4991 | 0.4982 | 0.4979 | 0.5153 | |
| Intel i9 10980xe | 18 cores | ||||||
| 192 GB Ram | 36 with hyper threading | ||||||
| Nvidia Quadro RTX 4000 | |||||||
| 8 GB ram | |||||||
It is example ch3305 from the book. We’ve updated the example since publication. A tar file is available on the fortranplus web site. Fortranplus - 4th edition *new* examples for interested people.
1 Like
I forgot to mention that the Cray run is on the Archer2 service. About ARCHER2
1 Like
asked an LLM to write a script to plot the data because I dislike tables, first by compiler comparison of the ops:
I
Then all ops per compiler
1 Like
I think it would be fair to add a note saying that ifort and ifx (that I know) translate do concurrent into OpenMP, and that one needs to activate that with flag. Otherwise it can be misleading and thus unfair.
And that also explains the gap for gfortran, as do concurrent so far is syntactically enabled but it hasn’t been mapped yet to OpenMP (AFAIK)
1 Like
It also makes clear that whole array operations are not necessarily (much) more expensive than explicit do-loops.
Hi hkvzjal the example comes from our openmp chapter, so all examples are compiled with openmp enabled.
Hi Arjen I think this helps support Damian’s statement that we should move from do loops to more modern constructs, as the cost difference isn’t very significant.
Thanks for the graphs. Looks good.
1 Like
And that’s exactly my point! You know it, but someone landing in this thread and just reading that result might not. I do believe that a little foot-note would be welcomed just for the sake of clarity. Such that people know that, if they don’t activate OpenMP they shall expect do concurrent to behave basically as plain old do, which happens with gfortran anyway because the compiler doesn’t do the translation yet.
And so, it brings me also to another question, in your nvidia results, I’m a bit surprised actually, because do concurrent is parallelized with OpenACC. Which flags did you use for the nvidia benchmarks? (… not surprised, actually the results are consistent, but it might need an extra row with do concurrent+OpenACC)
Exactly 
. I am a big fan of such constructs and I was surprised to learn of the almost transparant parallellisation that DO CONCURRENT offers.
Still, there is also a lot to learn about the gory details that surround GPUs. I have so far been unable to dig up publications that go beyond the mere syntax or one or two recommendations. (I have half a mind - or perhaps a quarter - to have a go at it myself.)
1 Like
Someone on HackNews commented the following (I do not necessarily agree):
The discussion there badly misunderstands the nature of ELEMENTAL procedures in Fortran and their relevance to parallel execution.
ELEMENTAL is relevant to DO CONCURRENT only indirectly. The ELEMENTAL attribute matters there only because it implies the PURE attribute by default, and PURE is required for procedures referenced in DO CONCURRENT. (Which is not a parallel construct, but that’s another matter.)
ELEMENTAL in array expressions (incl. FORALL) should not be understood as being a way for one procedure invocation to receive and return entire arrays as arguments and results. That would require buffering during the evaluation of an array expression. Instead, ELEMENTAL should be viewed (and implemented) as a means of allowing a function to be called as part of the implementation of unbuffered array expression execution.
ELEMENTAL has its roots in compilation for true vector machines. It once caused a function to have multiple versions generated: a normal one with scalar arguments, and a vector one with vector register arguments. This would allow a user-written ELEMENTAL function to be called in a vectorized DO loop, just like an intrinsic vectorizable function like SIN. A compiler for today’s SIMD “vector” ISAs could implement ELEMENTAL in a similar fashion.
1 Like