# Fortran is dead – Long live Fortran! **URL:** \ **Category:** Tutorials\ **Created:** [March 19, 2024, 9:27pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654 "2024-03-19T21:27:59Z")\ **Posts on this page:** 10\ **Page:** 1
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [March 19, 2024, 9:27pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/1 "2024-03-19T21:27:59Z")
[Torsten Hoefler](http://htor.ethz.ch/)’s random access spontaneous talk given at the 42nd anniversary Salishan Conference on High-Speed Computing in 2023. Discusses how to lift Fortran code to a data-centric representation to optimize it for accelerator devices. Work led by Alexandru Calotoiu in SPCL. [![](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/9/94aa8d62587fc71d903b1aad651cc46e2dee87af.jpeg "Fortran is dead – Long live Fortran!") ](https://www.youtube.com/watch?v=1s8TT3jE8Tg) ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [March 20, 2024, 12:28am UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/2 "2024-03-20T00:28:26Z")
For the loop nest at 8:05 Torsten says that not a single compiler can parallelize that loop automatically, but I remain confused if he’s talking about vectorization of the inner loop, or doing the reduction on the outer loop in parallel The loop appears to be taken out of the [dwarf-p-cloud](https://github.com/ecmwf-ifs/dwarf-p-cloudsc/tree/main) miniapp, specifically from lines 2655 - 2658 of the file [`cloudsc.F90`](https://github.com/ecmwf-ifs/dwarf-p-cloudsc/blob/main/src/cloudsc_fortran/cloudsc.F90) (which contains a single ~ 2800 line long subroutine 👀 ): ```fortran ! Backsubstitution ! step 1 DO JN=2,NCLV DO JM = 1,JN-1 ZQXN(KIDIA:KFDIA,JN)=ZQXN(KIDIA:KFDIA,JN)-ZQLHS(KIDIA:KFDIA,JN,JM) & & *ZQXN(KIDIA:KFDIA,JM) ENDDO ENDDO ``` Rewriting the array expression as a loop, ```fortran DO JM = 1,JN-1 DO JL = KIDIA, KFDIA ZQXN(JL,JN) = ZQXN(JL,JN) - factor*ZQXN(JL,JM) ENDDO ENDDO ``` gives a loop nest resembling what is shown in the slide. ---
**Author:** ![jorgeg](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jorgeg/32/6835_2.png) [@jorgeg](https://fortran-lang.discourse.group/u/jorgeg)\ **Post date:** [June 26, 2025, 10:22am UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/3 "2025-06-26T10:22:23Z")
I just saw this talk today and I wondered the same thing… ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [June 26, 2025, 4:34pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/4 "2025-06-26T16:34:01Z")
I think the problem for auto-vectorization is that the compiler might not be able to deduce there is no aliasing, because the `JM` terminates at `JN - 1`, and the results are accumulated into column `JN`. One way to fix this would be using an `IVDEP` directive, or using an `!$omp simd` directive. Edit: `do concurrent` (DC) would also work (assuming the Fortran compiler sees DC as an IVDEP annotation and not some other parallelization mode). ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [June 26, 2025, 4:42pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/5 "2025-06-26T16:42:08Z")
A very similar problem occurs during the trailing panel update in LU factorization. Levesque & Williamson write about this in “A Guidebook to Fortran on Supercomputers”: > In the next few examples we will examine cases where the compiler may not know if a loop is recursive. Loop nest 42030 is a traditionally coded Gaussian elimination scheme. Most compilers will try to vectorize the inner loop where row subtraction is taking place. Note that `A(J,K)`is a function of `A(J,I)`. Although both `I` and `J` are fixed values within that inner loop a compiler must attempt to determine whether `K` ever assumes the value of `I` (on other than the last iteration), since this would be recursive and must not be vectorized. > > There are two approaches for the compiler in optimizing this case. The first approach is to be smart enough to recognize that there is no possibility for recursion, since the initial (and lowest) value of `K` is `I+1`. The other approach is to generate one code sequence for executing the loop in scalar mode and another sequence to execute it in vector mode, then test at runtime whether `K` will assume the value of `I` in the inner loop. This second approach is called conditional compilation and is used quite extensively by CFT 1.15 [Cray Fortran]. > > To assure that the compiler does fully vectorizes the restructured loop 42031 we have inserted directives (for Alliant, Cray, and NEC in this case) that inform the compiler to ignore potential recursion within the loop. Figure 4.10 illustrates that the performance improves by 50-100% for longer vector lengths. The performance gain results from the elimination of the execution time test by the compiler. ```fortran C THE ORIGINAL C GAUSS ELIMINATION DO 43020 I = 1, MATDIM A(I,I) = 1. / A(I,I) DO 43020 J = I+1, MATDIM A(J,I) = A(J,I) * A(I,I) DO 43020 K = I+1, MATDIM A(J,K) = A(J,K) - A(J,I) * A(I,K) 43020 CONTINUE C THE RESTRUCTURED C GAUSS ELIMINATION DO 43021 I = 1, MATDIM A(I,I) = 1. / A(I,I) DO 43021 J = I+1, MATDIM A(J,I) = A(J,I) * A(I,I) CVD$ NODEPCHK CDIR$ IVDEP *VDIR NODEP DO 43021 K = I+1, MATDIM A(J,K) = A(J,K) - A(J,I) * A(I,K) 43021 CONTINUE ``` This corresponds to the BLAS routine [DGER](https://netlib.org/lapack/explore-html-3.6.1/d7/d15/group __double__ blas__level2_ga458222e01b4d348e9b52b9343d52f828.html#ga458222e01b4d348e9b52b9343d52f828) used by [DGETF2](https://netlib.org/lapack/explore-html-3.6.1/dd/d9a/group__double_g_ecomputational_ga8360b5b2c819e19c82bfd7e6b8285f74.html#ga8360b5b2c819e19c82bfd7e6b8285f74). ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [June 27, 2025, 1:09am UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/6 "2025-06-27T01:09:47Z")
> [@ivanpribec](#): > > […] we have inserted directives (for Alliant, Cray, and NEC in this case) that inform the compiler to ignore potential recursion within the loop. I wrote fortran code for all of those machines in the 1980s and early 1990s, and I used those same compiler directives many times. I think is this case, a slight rewrite of the code allows the compiler to vectorize the innermost loop. ```auto DO 43020 I = 1, MATDIM A(I,I) = 1. / A(I,I) DO 43020 J = I+1, MATDIM AJI = A(J,I) * A(I,I) A(J,I) = AJI DO 43020 K = I+1, MATDIM A(J,K) = A(J,K) - AJI * A(I,K) 43020 CONTINUE ``` The compiler will likely not even store the `AJI` variable in memory, it will just live in a register for its lifetime. Of course, one would like for a fortran compiler to recognize all cases of vectorization, or the lack of recursion in this case, but especially in the 1980s that did not always occur. However, if this is still the situation in the language almost 40 years later, perhaps there should be some standard way to tell the compiler what to do. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [June 27, 2025, 7:08am UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/7 "2025-06-27T07:08:48Z")
> [@RonShepard](#): > > However, if this is still the situation in the language almost 40 years later, perhaps there should be some standard way to tell the compiler what to do. `do concurrent`? Array expressions? `!$omp simd` also belongs to _a_ standard. If I remember correctly, the Intel Fortran compiler has `-qopenmp-simd` turned ON by default. The compiler would need to prove that `I+1 .. MATDIM` and `I` are disjoint sets, hence the loop can be safely vectorized. I guess there is a name for this type of reasoning (loop aliasing analysis?). Maybe @certik knows? * * * Here is a slight variation of Levesque’s loop and the DGER routine: ```fortran subroutine trailing_update(m,n,A,lda,k) integer, parameter :: dp = kind(1.0d0) integer, intent(in) :: m, n, lda, k real(dp), intent(inout) :: A(lda,n) integer :: i, j do j = k+1,n do i = k+1,m A(i,j) = A(i,j) - A(i,k)*A(k,j) end do end do end subroutine ``` When compiling with `gfortran -O2 -ftree-vectorize -march=skylake`, it creates both a scalar and a vector path (just like Cray CFT did): ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/b/b4b5a8b26f4aea7dff26c369f69c193913e24870.png) If you add `!GCC$ IVDEP` or `!$omp simd` the scalar path dissappears: ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/7/7bbf7d4f8c6e385e8450c44c7dac219d7f4c9d5a.png) Even the Intel compiler (`ifx -O3 -xSKYLAKE`) struggles to create a vector version. Instead it unrolls 4 iterations and uses scalar instructions: ```txt .LBB0_5: vmovsd xmm0, qword ptr [r13 + 8*rsi - 8] vmovsd xmm1, qword ptr [r11 + rbp - 24] vfnmadd213sd xmm1, xmm0, qword ptr [r14 + rbp - 24] vmovsd qword ptr [r14 + rbp - 24], xmm1 vmovsd xmm0, qword ptr [r13 + 8*rsi - 8] vmovsd xmm1, qword ptr [r11 + rbp - 16] vfnmadd213sd xmm1, xmm0, qword ptr [r14 + rbp - 16] vmovsd qword ptr [r14 + rbp - 16], xmm1 vmovsd xmm0, qword ptr [r13 + 8*rsi - 8] vmovsd xmm1, qword ptr [r11 + rbp - 8] vfnmadd213sd xmm1, xmm0, qword ptr [r14 + rbp - 8] vmovsd qword ptr [r14 + rbp - 8], xmm1 vmovsd xmm0, qword ptr [r13 + 8*rsi - 8] vmovsd xmm1, qword ptr [r11 + rbp] vfnmadd213sd xmm1, xmm0, qword ptr [r14 + rbp] vmovsd qword ptr [r14 + rbp], xmm1 add rbp, 32 cmp r15, rbp jne .LBB0_5 ``` When you add `!dir$ ivdep` or `!$omp simd`, then `ifx` also does the right thing: ```txt vbroadcastsd ymm1, xmm0 xor ebx, ebx .LBB0_7: vmovupd ymm2, ymmword ptr [r9 + 8*rbx] vfnmadd213pd ymm2, ymm1, ymmword ptr [r10 + 8*rbx] vmovupd ymmword ptr [r10 + 8*rbx], ymm2 add rbx, 4 cmp rbx, r8 jb .LBB0_7 ``` One more difference worthy of notice, ✅ the Intel compiler broadcasts the `A(k,j)` value before the hot loop . ❌ `gfortran` does the scalar-to-vector broadcast on each iteration (despite the fact that `rax` is constant in that loop). What is interesting is `IFORT` vectorizes this loop nest even without directives: ```auto ..B1.14: # Preds ..B1.12 ..B1.10 vbroadcastsd ymm3, xmm1 #10.38 vbroadcastsd ymm2, xmm0 #10.38 movsxd rcx, esi #9.9 mov r8, QWORD PTR [-48+rsp] #9.9[spill] ..B1.15: # Preds ..B1.15 ..B1.14 vmovupd ymm4, YMMWORD PTR [r8+rdx*8] #10.31 vfnmadd213pd ymm4, ymm3, YMMWORD PTR [rbx+rdx*8] #10.13 vmovupd YMMWORD PTR [rbx+rdx*8], ymm4 #10.13 vmovupd ymm5, YMMWORD PTR [r8+rdx*8] #10.31 vfnmadd213pd ymm5, ymm2, YMMWORD PTR [rdi+rdx*8] #10.13 vmovupd YMMWORD PTR [rdi+rdx*8], ymm5 #10.13 vmovupd ymm6, YMMWORD PTR [32+r8+rdx*8] #10.31 vfnmadd213pd ymm6, ymm3, YMMWORD PTR [32+rbx+rdx*8] #10.13 vmovupd YMMWORD PTR [32+rbx+rdx*8], ymm6 #10.13 vmovupd ymm7, YMMWORD PTR [32+r8+rdx*8] #10.31 vfnmadd213pd ymm7, ymm2, YMMWORD PTR [32+rdi+rdx*8] #10.13 vmovupd YMMWORD PTR [32+rdi+rdx*8], ymm7 #10.13 vmovupd ymm8, YMMWORD PTR [64+r8+rdx*8] #10.31 vfnmadd213pd ymm8, ymm3, YMMWORD PTR [64+rbx+rdx*8] #10.13 vmovupd YMMWORD PTR [64+rbx+rdx*8], ymm8 #10.13 vmovupd ymm9, YMMWORD PTR [64+r8+rdx*8] #10.31 vfnmadd213pd ymm9, ymm2, YMMWORD PTR [64+rdi+rdx*8] #10.13 vmovupd YMMWORD PTR [64+rdi+rdx*8], ymm9 #10.13 vmovupd ymm10, YMMWORD PTR [96+r8+rdx*8] #10.31 vfnmadd213pd ymm10, ymm3, YMMWORD PTR [96+rbx+rdx*8] #10.13 vmovupd YMMWORD PTR [96+rbx+rdx*8], ymm10 #10.13 vmovupd ymm11, YMMWORD PTR [96+r8+rdx*8] #10.31 vfnmadd213pd ymm11, ymm2, YMMWORD PTR [96+rdi+rdx*8] #10.13 vmovupd YMMWORD PTR [96+rdi+rdx*8], ymm11 #10.13 add rdx, 16 #9.9 cmp rdx, rcx #9.9 jb ..B1.15 # Prob 82% ``` Roughly speaking it interleaves the updates for two iterations in the second dimension (notice the `vbroadcastsd`), and unrolling 16 iterations along the first array dimension. It also uses different registers, presumably to avoid (false) instruction scheduling dependencies. I wonder what heuristics it uses to determine these unrolling factors. Obviously register availability is one. `ifx` does something similar as IFORT, when you use ```fortran do concurrent (j = k+1:n, i=k+1:m) A(i,j) = A(i,j) - A(i,k)*A(k,j) end do ``` But instead of chunks 16 x 2, it iterates across chunks of size 4 x 8. Also note the repeates use of `ymm0`: ```auto vmovapd xmmword ptr [rsp + 48], xmm8 vbroadcastsd ymm8, xmm8 vbroadcastsd ymm9, xmm1 vbroadcastsd ymm10, xmm2 vbroadcastsd ymm11, xmm3 vbroadcastsd ymm12, xmm4 vbroadcastsd ymm13, xmm5 vbroadcastsd ymm14, xmm6 vbroadcastsd ymm15, xmm7 xor eax, eax .LBB0_9: vmovupd ymm0, ymmword ptr [rbp + 8*rax] lea rdx, [r11 + 8*rax] vfnmadd213pd ymm0, ymm8, ymmword ptr [r11 + 8*rax] vmovupd ymmword ptr [r11 + 8*rax], ymm0 vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm9, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 add rdx, rcx vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm10, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 lea rdx, [rdx + rcx] vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm11, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 lea rdx, [rdx + rcx] vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm12, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 lea rdx, [rdx + rcx] vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm13, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 add rdx, rcx vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm14, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 add rdx, rcx vmovupd ymm0, ymmword ptr [rbp + 8*rax] vfnmadd213pd ymm0, ymm15, ymmword ptr [rcx + rdx] vmovupd ymmword ptr [rcx + rdx], ymm0 add rax, 4 cmp rax, r13 jbe .LBB0_9 ``` So three compilers, three different strategies. Two compilers (`gfortran` and `ifx`) only vectorize when the programmer asserts it is safe. Addendum: gfortran moves the broadcast before the loop if written like @RonShepard suggested: ```Fortran do j = k+1,n tmp = A(k,j) do i = k+1,m A(i,j) = A(i,j) - A(i,k)*tmp end do end do ``` ```txt vbroadcastsd ymm1, xmm2 xor edi, edi .L6: vmovupd ymm0, YMMWORD PTR [rcx+rdi] vfnmadd213pd ymm0, ymm1, YMMWORD PTR [rax+rdi] vmovupd YMMWORD PTR [rax+rdi], ymm0 add rdi, 32 cmp rdx, rdi jne .L6 ``` The iteration space is traversed in chunks of 4 x 1, same as `ifx` did. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [June 27, 2025, 3:09pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/9 "2025-06-27T15:09:01Z")
> [@ivanpribec](#): > > The compiler would need to prove that `I+1 .. MATDIM` and `I` are disjoint sets, hence the loop can be safely vectorized. This is obvious to a human, but I have no idea how optimizers deduce such things, or whether this is part of the language-specific front end or part of the general back end of a modern compiler. I did think about do concurrent, but my feeling is that construct is more about concurrent stripmining and subblocking that it is about recursion. However, the lack of recursion is indeed part of the semantics of that construct, so maybe that is sufficient? I also considered OpenMP, and I think that also works, but it implies multithreaded execution, which is also above and beyond just telling the compiler that no recursion is occurring within the loop. And, of course, OpenMP is not part of the language standard, it is a separate standard. Array syntax actually does the opposite of what the programmer needs; it tells the compiler if it can’t figure it out, then allocate temporary arrays for the operands and for the result and do everything with those temporary arrays. In this case, with just two floating point operations per loop cycle that workspace allocation and data movement would be a performance disaster. Thanks much for your effort within compiler explorer. I don’t understand the details of the machine instructions, but I think I can follow your discussion and understand the consequences of the various compiler choices. ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [June 27, 2025, 3:10pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/10 "2025-06-27T15:10:41Z")
> [@ivanpribec](#): > > The compiler would need to prove that `I+1 .. MATDIM` and `I` are disjoint sets, hence the loop can be safely vectorized. I guess there is a name for this type of reasoning (loop aliasing analysis?). Maybe @certik knows? I don’t know either yet. After beta we’ll focus on optimizations more with LFortran. I have many ideas though. 🙂 ---
**Author:** ![wspector](https://avatars.discourse-cdn.com/v4/letter/w/47e85d/32.png) [@wspector](https://fortran-lang.discourse.group/u/wspector)\ **Post date:** [June 27, 2025, 4:50pm UTC](https://fortran-lang.discourse.group/t/fortran-is-dead-long-live-fortran/7654/11 "2025-06-27T16:50:23Z")
> [@ivanpribec](#): > > There are two approaches for the compiler in optimizing this case. The first approach is to be smart enough to recognize that there is no possibility for recursion, since the initial (and lowest) value of `K` is `I+1`. The other approach is to generate one code sequence for executing the loop in scalar mode and another sequence to execute it in vector mode, then test at runtime whether `K` will assume the value of `I` in the inner loop. This second approach is called conditional compilation and is used quite extensively by CFT 1.15 [Cray Fortran]. IIRC, CFT 1.14 could also do it. At the time, Cray was also developing the next generation compiler - CFT77. There was a period where CFT could do tricks that CFT77 couldn’t - and vice versa. CFT77 was a globally optimizing compiler, whereas CFT mostly optimized within basic blocks, and towards the end, extended basic blocks. So CFT77 generated better scalar code from the start. But eventually CFT77 learned all of CFT’s vector tricks. What was really fun was when one of the compiler developers grafted the C compiler front end onto the CFT77 optimization/vectorization and code generation. It was one of those little “after midnight” projects which worked a little too well, and ended up causing a bit of a stir…