Как дать указание gfortran скомпилировать цикл, используя только векторизованные (максимально доступные SIMD) инструкции для распределяемых массивов?

#fortran #gfortran #simd #intel-fortran #avx2

# #fortran #gfortran #simd #intel-fortran #avx2

Вопрос:

Когда gfortran векторизирует уравнение, такое как x = y * z , компилятор будет использовать только векторные регистры (такие как YMM) и векторные коды операций, тогда и только тогда, когда компилятор знает, что все данные (x, y, z) могут быть загружены в указанные векторные регистры. Однако, если компилятор не знает, все ли данные могут быть загружены в регистры, он сгенерирует избыточный код, чтобы можно было оперировать оставшимися данными.

Этот избыточный код обычно не будет использовать максимально доступные векторные регистры SIMD и будет раздувать исполняемый файл. Я не проводил тщательного тестирования, но общие рассуждения наводят меня на мысль, что это может привести к пропускам ICACHE или предотвратить оптимизацию встраивания функций / подпрограмм.

Я хотел бы удалить этот избыточный код, поскольку я знаю, что все мои данные (x, y, z) будут идеально вписываться в векторные регистры YMM.

Ранее я более подробно описал свои наблюдения здесь: https://www.cfd-online.com/Forums/main/231759-fortran-assembly-how-remove-redundant-non-vectorized-code.html

Тем не менее, я хотел бы привести здесь простые примеры, чтобы продемонстрировать свою точку зрения.

Для операции: x(:) = y*z где x, y, z — все распределяемые массивы, компилятор не знает, можно ли полностью векторизовать операцию, поскольку он не знает, какова длина массивов. Если x, y, z являются целочисленными массивами и имеют длину, кратную 8, мы можем сделать вывод, что вся операция может быть выполнена просто с помощью векторизованных операций AVX2. Однако, поскольку длины массивов неизвестны во время компиляции, компилятору приходится генерировать избыточный код, который будет работать с массивами произвольной длины.

ССЫЛКА НА GODBOLT: https://gcc.godbolt.org/z/1fdzK8

Вот сборка, сгенерированная для x(:) = y*z части:

 .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
 

Мы это ясно видим .Раздел L7 приведенного выше кода содержит векторизованный код, а раздел .L6 содержит избыточный код, гарантирующий, что массивы произвольной длины будут работать с одним и тем же кодом.

Хотя это хорошая функция, я хотел бы надежно отключить ее и дать указание компилятору генерировать только векторизованный код в разделе .L7. Я не проводил тщательного тестирования производительности, однако я все равно хотел бы удалить раздел .L6.

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
 

Я также смог заставить компилятор gfortran генерировать только раздел векторизованного кода, используя некоторый компилятор voodoo. Мы просто используем do i=1,AND(N,-8) для указания компилятору, что цикл заканчивается на кратном 8. Это происходит, поскольку математически AND(N,-8) == N - MOD(N,8)

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

ССЫЛКА НА GODBOLT: https://gcc.godbolt.org/z/MW7cPG

Сгенерированная сборка была :

 .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
 

Однако это ненадежный метод для gfortran, и в некоторых моих тестовых примерах он не удаляет избыточное раздувание. AFAIK gfortran не имеет какой-либо директивы компилятора like !DIR$ ASSUME (MOD(N,8) .EQ. 0) , поэтому этот метод кажется мне крайне ненадежным.

Я хочу сделать эту оптимизацию надежным и последовательным способом.

ДОПОЛНИТЕЛЬНЫЕ ДАННЫЕ :

Те же операции были выполнены над, казалось бы, сложным кодом :

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

ССЫЛКА НА GODBOLT: https://gcc.godbolt.org/z/oxcWzd

В приведенном выше примере генерируемое раздувание довольно велико.

Вышеупомянутая оптимизация работает в компиляторе Fortran от Intel.

ССЫЛКА НА GODBOLT: https://gcc.godbolt.org/z/Y5orh1

Вышеупомянутая оптимизация работает в компиляторе gfortran.

ССЫЛКА НА GODBOLT: https://gcc.godbolt.org/z/1or59d

Комментарии:

1. Не могу помочь с проблемой векторизации, но отмечу, что если у вас есть выбор между a real и integer exponent , вероятно, хорошей идеей будет использовать an integer . Некоторые компиляторы могут не распознавать x**2.0 то же самое, x**2 что и, и, следовательно, он может вычислять мощность с помощью pow(x,2.0) or exp(2.0*log(x)) . Да, я знаю 2.0 , что это точно представимо.

2. Да, я провел сравнительный анализ между вещественным и целочисленным показателем. Целочисленные показатели значительно быстрее и хорошо векторизуются, в то время как реальные показатели создают скалярный код <<< ССЫЛКА >>>

3. Прочитав это несколько раз, я должен признать, что все еще думаю «почему?». Я бы провел некоторый сравнительный анализ, прежде чем вводить такое запутывание кода… Промахи в кэше инструкций вряд ли будут проблемой, поскольку, если код никогда не выполняется, он никогда не загружается в кэш (за исключением нечетного слова здесь и там, разделяющего строку кэша), и меня также не убеждают комментарии по встраиванию без каких-либо подтверждающих доказательств.

Ответ №1:

Fortran — по философии — не хочет, чтобы пользователь возился с низкоуровневой оптимизацией, потому что предполагается, что язык позволяет это. В вашем случае полная векторизация может быть достигнута путем определения объема как elemental функции:

   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
 

а затем выполнить цикл с использованием forall или 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
 

Это создает чистый, векторизованный код.

Обратите внимание, что do concurrent здесь включена векторизация процессора, но в некоторых будущих версиях может скрываться интересные вещи, такие как выгрузка на графический процессор.

Комментарии:

1. Для элементарной функции вы можете пропустить do цикл и просто указать vols = triVolume(x1,x2,x3,y1,y2,y3)

2. Графические процессоры уже здесь: developer.nvidia.com/blog /… , при условии , что вы можете заставить свой код работать с компилятором NVFORTRAN … Обратите внимание, что forall также устарел

3. @IanBush Я пропускаю любые упоминания о передаче памяти в и из GPU и пропускной способности.