Фортран: пишем параллельные программы03.04.2023 20:16

Современный Фортран представляет собой специализированный язык программирования, предназначенный в основном для написания вычислительных программ для векторно-конвейерных и параллельных архитектур. Эволюция стандартов языка Фортран была рассмотрена в предыдущих статьях — здесь и здесь.

На данный момент действующим стандартом языка Фортран является стандарт ISO 2018 года «Fortran 2018», готовится к принятию стандарт 2023 года. К сожалению, различные компиляторы Фортрана поддерживают требования стандартов в различной степени.

В этой статье мы попробуем написать простейшую параллелизуемую программу на языке Фортран, используя для этого методы конвейеризации и симметричной параллелизации и сравним их между собой, применив наиболее популярные компиляторы GNU Fortran и Intel Fortran.

В целом, компилятор Intel Fortran гораздо более полно реализует стандарт Fortran 2018. В частности, он поддерживает все имеющиеся в стандарте средства параллельных вычислений, в то время как GNU Fortran реализует только самые базовые из них (чего, впрочем, в ряде случаев более чем достаточно). С другой стороны, Intel Fortran, в отличие от GNU Fortran, не обеспечивает реализацию символьного типа CHARACTER (KIND=4) с поддержкой кодировки UCS-4, что может затруднить обработку не-ASCII текстов. Бытует мнение, что Intel Fortran обладает более мощным оптимизатором.

Постановка задачи

Напишем простейшую программу для реализации классического клеточного автомата игры «Жизнь». Не будем сейчас париться с вводом и выводом, исходную конфигурацию зададим в самой программе, а результирующую конфигурацию после заданного числа шагов выведем в файл. Нас будут интересовать сами вычислительные шаги клеточного автомата.

Для тестов, чтобы далеко не ходить, используется компьютер Mac mini с процессором Intel Core i3 с 4 физическими ядрами. Компиляторы GNU Fortran 12.2.0 и Intel Classic Fortran 2021.8.0 20221120.

Последовательная программа

Для начала напишем программу в чисто последовательном стиле. Напишем всё в одном файле, чтобы оптимизатору было легче работать.

program life
! чисто последовательный вариант программы

implicit none

! здесь мы задаём количество байтов
! для каждой ячейки - вдруг операции над
! целыми словами окажутся эффективными? (нет)
integer, parameter :: matrix_kind = 1

integer, parameter :: generations = 2 ! автомат рассматривает 2 поколения
integer, parameter :: rows = 1000, cols = 1000 ! размеры поля
integer, parameter :: steps = 10000 ! количество шагов

! описываем игровое поле. значения элементов могут быть целыми 0 или 1
integer (kind=matrix_kind) :: field (0:rows+1, 0:cols+1, generations)

integer :: thisstep = 1, nextstep =2 ! индексы массива для шагов
! при желании это можно легко обобщить на автомат с памятью больше 1 шага

integer :: i ! счётчик шагов
integer :: clock_cnt1, clock_cnt2, clock_rate ! для работы с таймером

! инициализируем поле на шаге thisstep начальной конфигурацией
call init_matrix (field (:, :, thisstep))

! засечём время
call system_clock (count=clock_cnt1)

! вызовем процедуру выполнения шага в цикле для заданного числа шагов
do i = 1, steps
  ! тут мы берём сечение массива по thisstep и преобразовываем в nextstep
  call process_step (field (:, :, thisstep), field (:, :, nextstep))
  ! следующий шаг становится текущим
  thisstep = nextstep
  ! а для следующего шага снова возвращаемся к другому сечению
  nextstep = 3 - thisstep
end do

! узнаем новое значение таймера и его частоту
call system_clock (count=clock_cnt2, count_rate=clock_rate)

! напечатаем затраченное время и оценку производительности
print *, (clock_cnt2-clock_cnt1)/clock_rate, 'сек, ', & 
  rows*cols*steps/(clock_cnt2-clock_cnt1)*clock_rate, 'ячеек/с'

! выведем результирующую конфигурацию в файл для контроля
call output_matrix (field (:, :, thisstep))

! разместим подпрограммы тут же, чтобы оптимизатору было проще
contains

! проинициализируем, просто воткнув одну "мигалку" в чистое поле
pure subroutine init_matrix (m)
  integer (kind=matrix_kind), intent (out) :: m (:,:)
  m = 0
  m (50, 50) = 1
  m (50, 51) = 1
  m (50, 52) = 1
end subroutine init_matrix

! выведем матрицу в файл при помощи пробелов, звёздочек и грязного хака
subroutine output_matrix (m)
  integer (kind=matrix_kind), intent (in) :: m (0:,0:)
  integer :: rows, cols
  integer :: i, j
  integer :: outfile
  rows = size (m, dim=1) - 2
  cols = size (m, dim=2) - 2
  open (file = 'life.txt', newunit=outfile)
  do i = 1, rows
    ! выводим в каждой позиции строки символ, код которого является
    ! суммой кода пробела и значения ячейки (0 или 1), умноженного
    ! на разность между звёздочкой и пробелом
    write (outfile, '(*(A1))') (char (ichar (' ') + &
      m(i, j)*(ichar ('*') - ichar (' '))), j=1, cols)
  end do
  close (outfile)
end subroutine output_matrix

! здесь самое интересное – обработка шага
! для начала простой последовательный алгоритм

pure subroutine process_step (m1, m2)

  integer (kind=matrix_kind), intent (in) :: m1 (0:,0:)
  integer (kind=matrix_kind), intent (out) :: m2 (0:,0:)
  integer :: rows, cols
  integer :: i, j, s

  ! восстанавливаем значения rows и cols
  ! конечно, мы могли бы из просто передать в параметрах, но так культурнее
  rows = size (m1, dim=1) - 2
  cols = size (m1, dim=2) - 2

  ! обычные последовательные вложенные циклы
  ! поскольку в Фортране массивы хранятся по столбцам, то j раньше i
  do j = 1, cols
    do i = 1, rows

      ! считаем количество живых соседей
      s = m1 (i-1, j) + m1 (i+1, j) + m1 (i-1, j-1) + m1 (i+1, j-1) + &
          m1 (i, j-1) + m1 (i-1, j+1) + m1 (i, j+1) + m1 (i+1, j+1)

      ! присваиваем значение выходной клетке 
      select case (s)
        case (3)
          m2 (i, j) = 1
        case (2)
          m2 (i, j) = m1 (i, j)
        case default
          m2 (i, j) = 0
      end select

    end do
  end do

  ! закольцуем игровое поле, используя гало в массиве, 
  ! дублирующее крайние элементы с другой стороны массива
  m2 (0,:) = m2 (rows, :)
  m2 (rows+1, :) = m2 (1, :)
  m2 (:, 0) = m2 (:, cols)
  m2 (:, cols+1) = m2 (:, 1)

end subroutine process_step

end program life

Откомпилируем нашу программу при помощи GNU Fortran и Intel Fortran:

$ gfortran life_seq.f90 -o life_seq_g -O3 -ftree-vectorize -fopt-info-vec -flto $ ifort life_seq.f90 -o life_seq -O3

Запустим:

$ ./life_seq_g 11 сек, 125172000 ячеек/с $ ./life_seq 14 сек, 94060000 ячеек/с

125 лямов в секунду у GNU Fortran против 94 лямов у Intel Fortran.

Сразу забежим вперёд и скажем, что Intel Fortran, вопреки своей репутации, отстаёт от GNU Fortran на всех дальнейших сопоставимых тестах. Возможно, ситуация была бы другой на более мощных и современных процессорах, а также она другая на массивно-параллельных системах, на причинах чего мы остановимся в следующей статье. Но пока вот так.

Давайте, может, попробуем 32-разрядные целые вместо байтов?

integer, parameter :: matrix_kind = 4

$ ./life_seq_g
13 сек, 107441000 ячеек/с
$ ./life_seq
19 сек, 70660000 ячеек/с

Как видим, ничего хорошего нам это не дало, а Intel Fortran ещё более подкосило.

1. Матричная программа

Некоторые люди думают, что, если заменить циклы неявными вычислениями с матрицами, то это невероятно оптимизирует код. Посмотрим, так ли это. Поменяем нашу любимую подпрограмму process_step:

! обработка шага операциями с матрицами

pure subroutine process_step (m1, m2)
  
  integer (kind=matrix_kind), intent (in) :: m1 (0:,0:)
  integer (kind=matrix_kind), intent (out) :: m2 (0:,0:)
  integer :: rows, cols
  integer s (0:size(m1,dim=1)-1, 0:size (m1,dim=2))
    
  rows = size (m1, dim=1) - 2
  cols = size (m1, dim=2) - 2

  ! вычислим матрицу s, которая повторяет по форме и размерам матрицу m1
  ! и содержит в каждом элементе количество живых соседей клетки

  s = m1(0:rows-1,:) + m1(2:rows+1,:) + m1(0:rows-1,0:cols-1) + & 
      m1(2:rows+1,0:cols-1) + m1(:,0:cols-1) + m1(0:rows-1,2:cols+1) + &          
      m1(:,2:cols+1) + m1(2:rows+1,2:cols+1)

  ! завернём края ещё до вычислений

  s (0,:) = s (rows, :)
  s (rows+1, :) = s (1, :)
  s (:, 0) = s (:, cols)  
  s (:, cols+1) = s (:, 1)

  ! и применим оператор матричной обработки where

  where (s==3 .or. s==2 .and. m1 == 1)
    m2 = 1
  elsewhere
    m2 = 0
  end where

end subroutine process_step

Вернёмся к matrix_kind = 1 и проверим мощь матричных операторов:

$ ./life_mat_g 11 сек, 126815000 ячеек/с $ ./life_mat 25 сек, 55580000 ячеек/с

Как видим, у GNU Fortran результат чуть-чуть хуже чисто последовательного алгоритма, а Intel Fortran вообще расстроился от таких дел.

При этом надо ещё отметить, что Intel Fortran по умолчанию размещает очень мало памяти для стека, и увеличение размеров игрового поля (а вместе с ним и размещаемой на стеке переменной s в матричном варианте) приводит к выпадению программы в кору. GNU Fortran свободно работает при настройках по умолчанию с огромным размером поля.

С другой стороны, складывается впечатление, что здесь можно серьёзно соптимизировать матричный алгоритм, чтобы не перебирать одни и те же элементы массива трижды при движении по матрице. Возможно, кто-то из читателей предложит своё решение.

2. SMP параллелизм через OpenMP

Обе предыдущие программы были чисто последовательными, хотя компиляторы немножко векторизовали операции. Это неинтересно. Давайте извлечём пользу из наличия нескольких ядер в процессоре, причём сделаем это самым простым и грубым способом — через OpenMP:

! обратите внимание, что подпрограмма, управляющая внутри себя 
! параллелизмом с помощью директив omp, не может быть объявлена чистой,
! так как это очевидный побочный эффект. декларация pure привела бы
! к ошибке компиляции

impure subroutine process_step (m1, m2)

  integer (kind=matrix_kind), intent (in) :: m1 (0:,0:)
  integer (kind=matrix_kind), intent (out) :: m2 (0:,0:)
  integer :: rows, cols
  integer :: i, j, s

  rows = size (m1, dim=1) - 2
  cols = size (m1, dim=2) - 2

  ! внешний цикл исполняется параллельно на ядрах SMP.
  ! переменные i и s свои в каждой параллельной ветке кода
  
  !$omp parallel do private (i, s)
  do j = 1, cols
    do i = 1, rows
      s = m1 (i-1, j) + m1 (i+1, j) + m1 (i-1, j-1) + m1 (i+1, j-1) + &
          m1 (i, j-1) + m1 (i-1, j+1) + m1 (i, j+1) + m1 (i+1, j+1)
      select case (s)
        case (3)
          m2 (i, j) = 1
        case (2)
          m2 (i, j) = m1 (i, j)
        case default
          m2 (i, j) = 0
      end select
    end do
  end do
 !$end parallel do
    
  m2 (0,:) = m2 (rows, :)
  m2 (rows+1, :) = m2 (1, :)
  m2 (:, 0) = m2 (:, cols)
  m2 (:, cols+1) = m2 (:, 1)
          
end subroutine process_step

Не забудем подключить OpenMP при компиляции:

$ gfortran life_omp.f90 -o life_omp_g -O3 -ftree-vectorize -fopt-info-vec -flto -fopenmp $ ifort life_omp.f90 -o life_omp -O3 -qopenmp

И запустим:

$ ./life_omp_g 3 сек, 377022000 ячеек/с $ ./life_omp 3 сек, 356690000 ячеек/с

Теперь наш цикл выполняется одновременно на 4 ядрах процессора, за счёт чего выполнение ускорилось в 3 с лишним раза. По-прежнему, однако, GNU Fortran чуть впереди Intel Fortran’а.

3. SMP параллелизм через DO CONCURRENT

Попробуем переписать нашу программу стандартными средствами параллельного SMP программирования языка Фортран, без использования внешнего API OpenMP:

! подпрограмма снова может быть чистой, так как она не управляет нитками

pure subroutine process_step (m1, m2)
  
  integer (kind=matrix_kind), intent (in) :: m1 (0:,0:)
  integer (kind=matrix_kind), intent (out) :: m2 (0:,0:)
  integer :: rows, cols
  integer :: i, j, s
  
  rows = size (m1, dim=1) - 2
  cols = size (m1, dim=2) - 2
  
  ! так выглядит параллельный цикл в стандарте Фортрана
  ! как и в OpenMP,здесь распараллелен только внешний цикл

  do concurrent (j = 1:cols) local (i, s)
    do i = 1, rows

      s = m1 (i-1, j) + m1 (i+1, j) + m1 (i-1, j-1) + m1 (i+1, j-1) + &
          m1 (i, j-1) + m1 (i-1, j+1) + m1 (i, j+1)+ m1 (i+1, j+1)
 
      select case (s)
        case (3)
          m2 (i, j) = 1
       case (2)
          m2 (i, j) = m1 (i, j)
        case default   
          m2 (i, j) = 0
      end select

    end do
  end do  
  
  m2 (0,:) = m2 (rows, :)
  m2 (rows+1, :) = m2 (1, :)
  m2 (:, 0) = m2 (:, cols)
  m2 (:, cols+1) = m2 (:, 1)

end subroutine process_step

Здесь нас ждёт некоторое разочарование, потому что конструкция DO CONCURRENT в GNU Fortran реализована мало и плохо. Предложение LOCAL не может быть оттранслировано этим компилятором. И даже если бы мы как-то вывернулись из этого положения, то GNU Fortran всё равно преобразует DO CONCURRENT в обычный последовательный цикл DO (в интернете встречаются утверждения, что иногда GNU Fortran способен распараллелить DO CONCURRENT, но автору не удалось достичь такого эффекта).

Поэтому трансляцию этого примера мы можем выполнить только в Intel Fortran (обратите внимание, что компилятору всё равно нужна многонитевая библиотека OpenMP для параллелизации, без неё цикл будет откомпилирован в последовательный код):

$ ifort life_con2.f90 -o life_con -O3 -qopenmp

Запустим:

$ ./life_con 3 сек, 355890000 ячеек/с

Этот результат лучше всего, что мы видели в Intel Fortran, хотя немного не дотягивает до результата GNU Fortran с OpenMP.

4. Больше SMP параллелизма

Синтаксис оператора DO CONCURRENT как бы намекает, что мы можем объединить внутренний и внешний циклы в один параллельный цикл по двум параметрам. Посмотрим, что это даст:

! подпрограмма снова может быть чистой, так как она не управляет нитками 
! объединяем циклы do в общий do concurrent

pure subroutine process_step (m1, m2)

  integer (kind=matrix_kind), intent (in) :: m1 (0:,0:)
  integer (kind=matrix_kind), intent (out) :: m2 (0:,0:)
  integer :: rows, cols
  integer :: i, j, s
  
  rows = size (m1, dim=1) - 2
  cols = size (m1, dim=2) - 2
  
  ! так выглядит параллельный цикл в стандарте Фортрана
  ! здесь распараллелен как внешний, так и внутренний цикл
  ! в единую параллельную конструкцию, параметризованную по j и i
  do concurrent (j = 1:cols, i = 1:rows) local (s)

    s = m1 (i-1, j) + m1 (i+1, j) + m1 (i-1, j-1) + m1 (i+1, j-1) + &
        m1 (i, j-1) + m1 (i-1, j+1) + m1 (i, j+1) + m1 (i+1, j+1)

    select case (s)
      case (3)
        m2 (i, j) = 1
      case (2)
        m2 (i, j) = m1 (i, j)
      case default
        m2 (i, j) = 0
    end select  

  end do
  
  m2 (0,:) = m2 (rows, :)
  m2 (rows+1, :) = m2 (1, :)
  m2 (:, 0) = m2 (:, cols)
  m2 (:, cols+1) = m2 (:, 1)

end subroutine process_step

Что же это нам даёт?

$ ./life_con2 4 сек, 308920000 ячеек/с

Компилятор увлёкся обилием возможностей и ухудшил результат. Так что параллелить всё же надо с умом.

Вывод

Мы рассмотрели компиляцию простейшей программы на современном Фортране с использованием средств векторизации и симметричных параллельных вычислений. В результате тестов Intel Fortran показал некоторое преимущество в поддержке возможностей языка, а GNU Fortran — в скорости работы кода. При этом, однако, не надо забывать, что Intel Fortran поддерживает мощные методы совместной оптимизации раздельно расположенных в исходных файлах единиц компиляции, поэтому для большой программы сравнительный результат мог бы быть другим.

Вероятно, в следующей статье мы постараемся рассмотреть средства поддержки массивно-параллельных архитектур, имеющиеся в современном Фортране.