How to instruct gfortran to compile a loop using only vectorized(highest available SIMD) instructions for allocatable arrays?

Question

When gfortran vectorizes an equation such as x = y*z , the compiler will only use vector registers (such as YMM) and vector opcodes, if and only if the compiler knows that all of the data (x,y,z) can be loaded on said vector registers. However if the compiler doesn't know whether all of the data can be loaded onto the registers, it will generate redundant code such that the remaining data can be operated upon.

This redundant code will generally not use the highest available SIMD vector registers and will bloat the executable. I haven't done extensive testing, but common reasoning leads me to believe that it could cause ICACHE misses or prevent function/subroutine inlining optimizations.

I would like to remove this redundant code since I know that all of my data (x,y,z) will fit perfectly on the YMM vector registers.

I have previously described my observations in greater details here : https://www.cfd-online.com/Forums/main/231759-fortran-assembly-how-remove-redundant-non-vectorized-code.html

However, I would like to give simple examples here to showcase my point.

For the operation : x(:) = y*z where x,y,z are all allocatable arrays, the compiler doesn't know if the operation can be completely vectorized since it doesn't know what is the length of the arrays. If x,y,z are integer arrays and have a length that is a multiple of 8, we can deduce that the whole operation can be done simply using AVX2 vectorized operations. However since the array lengths aren't known at compile time, the compiler has to generate redundant code that will work on arrays of arbitrary lengths.

GODBOLT LINK : https://gcc.godbolt.org/z/1fdzK8

Here's the assembly generated for the x(:) = y*z portion:

.L7:
        vmovdqu ymm1, YMMWORD PTR [r13+0+rax]
        vpmulld ymm0, ymm1, YMMWORD PTR [r12+rax]
        vmovdqu YMMWORD PTR [r14+rax], ymm0
        add     rax, 32
        cmp     rdx, rax
        jne     .L7
        mov     rdx, r15
        and     rdx, -8
        lea     rax, [rdx+1]
        cmp     rdx, r15
        je      .L22
        vzeroupper
.L6:
        lea     rdx, [rax-1]
        mov     esi, DWORD PTR [r13+0+rdx*4]
        imul    esi, DWORD PTR [r12+rdx*4]
        mov     DWORD PTR [r14+rdx*4], esi
        lea     rdx, [rax+1]
        cmp     r15, rdx
        jl      .L10
        mov     esi, DWORD PTR [r13+0+rax*4]
        imul    esi, DWORD PTR [r12+rax*4]
        mov     DWORD PTR [r14+rax*4], esi
        lea     rsi, [rax+2]
        cmp     r15, rsi
        jl      .L10
        mov     edi, DWORD PTR [r13+0+rdx*4]
        imul    edi, DWORD PTR [r12+rdx*4]
        mov     DWORD PTR [r14+rdx*4], edi
        lea     rdx, [rax+3]
        cmp     r15, rdx
        jl      .L10
        mov     edi, DWORD PTR [r13+0+rsi*4]
        imul    edi, DWORD PTR [r12+rsi*4]
        mov     DWORD PTR [r14+rsi*4], edi
        lea     rsi, [rax+4]
        cmp     r15, rsi
        jl      .L10
        mov     edi, DWORD PTR [r13+0+rdx*4]
        imul    edi, DWORD PTR [r12+rdx*4]
        mov     DWORD PTR [r14+rdx*4], edi
        lea     rdx, [rax+5]
        cmp     r15, rdx
        jl      .L10
        mov     edi, DWORD PTR [r13+0+rsi*4]
        add     rax, 6
        imul    edi, DWORD PTR [r12+rsi*4]
        mov     DWORD PTR [r14+rsi*4], edi
        cmp     r15, rax
        jl      .L10
        mov     eax, DWORD PTR [r12+rdx*4]
        imul    eax, DWORD PTR [r13+0+rdx*4]
        mov     DWORD PTR [r14+rdx*4], eax

We can clearly see that .L7 section of the code above contains the vectorized code and the .L6 section contains redundant code to ensure that arrays of arbitrary lengths will work with the same code.

While this is a good feature, I would like to reliably turn it off and instruct the compiler to only generate the vectorized code in the .L7 section. I haven't done extensive performance testing, however I would still like to remove the .L6 section.

I was able to make Intel's Fortran compiler to only generate the vectorized code section using compiler directives such as:

!DIR$ ASSUME_ALIGNED X: 32
!DIR$ ASSUME (MOD(N,8) .EQ. 0)
do i=1,n
    x(i) = y(i)*z(i)
enddo

GODBOLT LINK : https://gcc.godbolt.org/z/PPf1zE

Generated assembly was :

..B1.15:                        # Preds ..B1.15 ..B1.14
        vmovdqu   ymm0, YMMWORD PTR [r13+r12*4]                 #22.12
        vpmulld   ymm1, ymm0, YMMWORD PTR [r14+r12*4]           #22.5
        vmovdqu   YMMWORD PTR [rax+r12*4], ymm1                 #22.5
        add       r12, 8                                        #21.1
        cmp       r12, rcx                                      #21.1
        jb        ..B1.15       # Prob 82%                      #21.1

I was also able to make gfortran compiler to only generate the vectorized code section by using some compiler voodoo. We simply use do i=1,AND(N,-8) to instruct to the compiler that the loop ends at a multiple of 8. It happens since mathematically AND(N,-8) == N - MOD(N,8)

do i=1,AND(n,-8)
    x(i) = y(i)*z(i)
end do

GODBOLT LINK : https://gcc.godbolt.org/z/MW7cPG

Generated assembly was :

.L19:
        mov     rcx, QWORD PTR [rsp+16]
        vmovdqu ymm4, YMMWORD PTR [rcx+rax]
        mov     rcx, QWORD PTR [rsp+8]
        vpmulld ymm0, ymm4, YMMWORD PTR [r15+rax]
        vmovdqu YMMWORD PTR [rcx+rax], ymm0
        add     rax, 32
        cmp     rax, rdx
        jne     .L19

However, this is not a reliable method for gfortran and fails to remove the redundant bloat in some of my test cases. AFAIK gfortran doesn't have any compiler directive like !DIR$ ASSUME (MOD(N,8) .EQ. 0), so this method seems highly unreliable to me.

I'm looking to do this optimization in a reliable and consistent way.

EXTRA DATA :

Same operations were performed on the seemingly complex code :

vols(:) = 0.5*sqrt((x1*(y2-y3))**2.0 + (x2*(y3-y1))**2.0 + (x3*(y1-y2))**2.0)

GODBOLT LINK : https://gcc.godbolt.org/z/oxcWzd

The bloat generated is quite high in the example above.

Aforementioned optimization works in Intel's Fortran compiler.

GODBOLT LINK : https://gcc.godbolt.org/z/Y5orh1

Aforementioned optimization works in gfortran compiler.

GODBOLT LINK : https://gcc.godbolt.org/z/1or59d

user_4702870 · Answer

Fortran - by philosophy - does not wish the user to tinker with low-level optimization, because the language is supposed to enable that. In your case, full vectorization can be achieved by defining volume as an elemental function:

  elemental real function triVolume(x1,x2,x3,y1,y2,y3) result(vol)
     real, intent(in) :: x1,x2,x3,y1,y2,y3
     vol = 0.5*sqrt((x1*(y2-y3))**2 + (x2*(y3-y1))**2 + (x3*(y1-y2))**2)
  end function triVolume

and then looping through it using forall or do concurrent:

! Calculate volume of a triangle using cross product
do concurrent (i=1:n)
    vols(i) = triVolume(x1(i),x2(i),x3(i),y1(i),y2(i),y3(i))
end do

That produces clean, vectorized code.

Note that do concurrent is here enabling CPU vectorization, but may, in some future releases, hide interesting stuff like off-loading to a GPU.

How to instruct gfortran to compile a loop using only vectorized(highest available SIMD) instructions for allocatable arrays?

1 Answers