В первой части статьи мы рассмотрели написание на современном Фортране простой программы, реализующей клеточный автомат "Жизнь", в виде классического последовательного кода (SISD), матричных операций (SIMD) и параллельных конструкций SMP (SIMD с частью функций MIMD). Сейчас мы будем рассматривать использование конструкций Фортрана для программирования массивно-параллельных архитектур (MPP), к которым, в частности, относятся современные суперкомпьютеры. Такие архитектуры реализуют классическую схему MIMD.

Как и в прошлой статье, мы используем для иллюстрации материала компиляторы GNU Fortran и Intel Fortran.

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

Как и в первой части статьи, мы продолжим реализовывать тот же самый клеточный автомат с теми же самыми входными и выходными данными. К сожалению, компилятор Intel Fortran не поддерживает используемую для программирования MPP архитектуру Coarray под операционной системой macOS, поэтому примеры для Intel Fortran мы в этот раз будем выполнять на другой машине с Linux. Это делает невозможным сравнение результатов в абсолютных цифрах, но по относительным оценкам всё по-прежнему будет ясно.

Для GNU Fortran, как и в первой части статьи, будет ипользоваться компилятор версии 12.2.0 и библиотека OpenCoarrays на компьютере Mac mini с 4-ядерным процессором Intel Core i3 @ 3.6 GHz под управлением macOS. Для Intel Fortran будут использоваться компилятор ifort версии  2021.9.0 20230302 и компилятор ifx 2023.1.0 20230320 на компьютере IBM x3250 M4 с 8-ядерным процессором Xeon 4C E3-1270v2 @ 3.9 GHz под управлением Linux. Ядер в этом процессоре больше, чем в Маке, но он постарее и потормознее.

0-4. Повторение темы

Напомним, что в прошлый раз мы достигли следующих результатов на Маке.

Последовательная программа в режиме автопараллелизации, gfortran и ifort:

$ ./life_seq_g

          11 сек,    124773000 ячеек/с

$ ./life_seq

           4 сек,    338120000 ячеек/с

Параллельная SMP программа средствами OpenMP, лучший результат gfortran:

$ ./life_omp_g

3 сек, 377022000 ячеек/с

$ ./life_omp

3 сек, 356690000 ячеек/с

Параллельная SMP программа средствами DO CONCURRENT, лучший результат ifort:

$ ./life_con

3 сек, 355890000 ячеек/с

Повторим эти результаты на нашем сервере Linux, чтобы можно было в дальнейшем опираться на какие-то сопоставимые значения (с автопараллелизацией ifort; ifx её не поддерживает):

$ ./life_seq_ifort

          9 сек,   144360000 ячеек/с

$ ./life_seq_ifx

          19 сек,     71030000 ячеек/с

$ ./life_con_ifort

           8 сек,  162910000 ячеек/с

$ ./life_con_ifx

          37 сек,     37730000 ячеек/с

Это в целом совпадает с результатами, приведёнными @mobi в комментариях к первой части статьи.

Прежде чем перейти к изложению нового материала, вспомним логику работы нашей программы в архитектуре SMP.

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

С другой стороны, архитектура SMP предопределила используемую нами модель оперативной памяти. Все объекты в памяти разделялись между всеми параллельными процессами, одинаково принадлежа им всем, за исключением некоторых отдельных локальных для процесса переменных, которые мы специально описывали как LOCAL или OMP PRIVATE.

5. MPP параллелизм через комассивы

Модель массивно-параллельных вычислений совершенно другая. Она просто излагается, но не всегда легко укладывается в голове для практического применения, поэтому тут нужно некоторое внимание.

Для начала, давайте разберёмся, чем вообще плохо SMP? Не особо парясь, парой подсказок компилятору мы можем ускорять свои последовательные программы в несколько раз – казалось бы, хватило бы только ядер, как пирату в мультфильме.

К сожалению, пулемётную ленту из ядер мы организовать не можем, так как очень быстро узким местом становится оперативная память. Разделение доступа к памяти между несколькими арифметико-логическими устройствами сложно чисто технически, и даже если мы это всё нужным образом припаяем и электрически нагрузим, то всё равно наши ядра будут конкурировать за доступ к одним и тем же страницам памяти (она же разделяемая) на более высоком уровне. Поэтому больше десятка-другого ядер в архитектуре SMP эффективно не подключить.

Что же нам делать, если мы хотим построить суперкомпьютер, включающий десятки тысяч вычислительных ядер (то есть собственно массивно-параллельную архитектуру)? В таком случае используются вычислительные узлы, каждый из которых содержит собственную оперативную память. Узлы объединены специальной высокоскоростной вычислительной сетью, позволяющей синхронизировать отдельные области памяти при помощи сообщений между узлами. Сам узел при этом чаще всего внутри себя, то есть на более низком уровне иерархии, использует несколько ядер, объединённых архитектурой SMP. Архитектура MPP не обязательно ограничивается двумя уровнями иерархии, как тут упрощённо описано; уровней может быть и больше (например, стойки суперкомпьютера взаимодействуют между собой с меньшей скоростью, чем узлы внутри стойки).

Программная архитектура MPP систем восходит к транспьютерам и ранним разработкам массивно-параллельного кода на языке Оккам.

Среда выполнения организована таким образом, что программа запускается одновременно в нескольких экземплярах, по одному экземпляру на узел (или ядро узла, что пока что для нас логически неважно). Текст и машинный код программы при этом на каждом используемом узле один и тот же, но динамика выполнения этого кода разная.

Так как код одинаков, то получается, что динамика его выполнения должна зависеть от каких-то внешних обстоятельств. И самым главным из таких обстоятельств является просто-напросто номер узла. Этот номер узла передаётся программе системной утилитой, запускающей массивно-параллельный код на массиве узлов. В Фортране номер узла, на котором в данный момент выполняется код, выдаётся функцией THIS_IMAGE() и всегда принимает значения от 1 до общего количества выполняющих программу узлов, которое, в свою очередь, выдаётся функцией NUM_IMAGES().

Как ясно из вышесказанного, каждый параллельно выполняющийся процесс в архитектуре MPP имеет свою собственную оперативную память. По этой причине все переменные, описанные в фортрановской программе, хранят собственные значения для каждого номера узла. В качестве исключения из этого правила, в коде могут быть специальным образом описаны комассивы, позволяющие программе получить доступ к оперативной памяти другого узла. Такой доступ к элементам комассивов имеет специальный синтаксис с квадратными скобками и не может применяться везде, где допустимо использование переменных в локальной оперативной памяти узла. Например, по очевидным причинам нельзя создать указатель на объект в чужой оперативной памяти.

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

Наш общий подход будет типичным для этого стиля программирования. Игровое поле мы объявим комассивом, к которому потенциально имеют доступ все узлы. Но фактически каждый узел будет заниматься только своей частью поля. Для этого разделим все столбцы матрицы field на полосы по числу узлов. Каждый узел рассчитывает значения элементов массива в своей полосе и заполняет полосу массива в своей оперативной памяти. Кроме того, по краям полосы узла находятся две единичные полоски, интересные его соседям. Значения в этих полосках узел копирует в оперативную память своих соседей.

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

Например, если у нас массив имеет 1002 столбца (1000 основных и 2 по краям для свёртки), а программа исполняется на 10 узлах, то узел номер 1 вычисляет столбцы 0:100, узел номер 2 – 101:200 и т.д. При этом столбец 100, рассчитанный узлом 1, помещается им не только в свою оперативную память, но и в память узла 2, так как этот столбец понадобится узлу 2 для расчёта столбца 101. Такие перекрывающиеся части комассивов называются “гало”.

Итак, вот новый код:

program life_mpp

    implicit none

    integer, parameter :: matrix_kind = 4
    integer, parameter :: generations = 2
    integer, parameter :: rows = 1000, cols = 1000
    integer, parameter :: steps = 10000

    integer (kind=matrix_kind) :: field (0:rows+1, 0:cols+1, generations) [*]
    integer :: thisstep = 1, nextstep =2
    integer :: i
    integer :: clock_cnt1, clock_cnt2, clock_rate

    integer, allocatable :: cols_lo (:), cols_hi (:) ! диапазоны столбцов для узлов

    integer :: me ! номер текущего узла, чтобы обращаться покороче

    me = this_image()

    print *, "it's me:", me

    ! заполним таблицы верхних и нижних границ полос для узлов

    allocate (cols_lo (num_images()), cols_hi (0:num_images()))
    cols_hi (0) = 0
    do i = 1, num_images()
      cols_lo (i) = cols_hi (i-1) + 1
      cols_hi (i) = cols * i / num_images()
    end do

    ! проинициализируем матрицу (инициализация тоже изменилась)

    call init_matrix (field (:, :, thisstep))

    sync all

    ! первый узел займётся хронометрированием

    if (me == 1) call system_clock (count=clock_cnt1)

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

    do i = 1, steps
      call process_step (field (:, :, thisstep), field (:, :, nextstep))
      thisstep = nextstep
      nextstep = 3 - thisstep
      sync all
    end do

    ! первый узел заканчивает хронометрирование, печатает результат,
    ! собирает матрицу и пишет в файл

    if (me == 1) then

      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, 'ячеек/с'

      ! мы хотим вывести матрицу в файл, а разные её столбцы хранятся
      ! на разных узлах. надо собрать всю матрицу на пишущем узле

      do i = 2, num_images()
        field (:, cols_lo(i):cols_hi(i), thisstep) = &
          field (:, cols_lo(i):cols_hi(i), thisstep) [i]
      end do

      call output_matrix (field (:, :, thisstep))

    end if

    ! надо синхронизироваться в конце, а то освободится память
    ! завершившегося процесса, пока её не скопировали
  
    sync all ! уходим всей бригадой

    contains

    pure subroutine init_matrix (m)

      integer (kind=matrix_kind), intent (out) :: m (0:,0:)
      integer j
 
      ! обнулим поле в своей полосе и гало
      ! не лезем далеко в чужие полосы, чтобы не распределять 
      ! ненужную узлу память

      do j = cols_lo(me)-1, cols_hi(me)+1
        m (:, j) = 0
      end do

      ! первый и последний узлы обнуляют кромки поля

      if (me == 1)             m (:, 0) = 0
      if (me == num_images())  m (:, cols+1) = 0

      ! нарисуем "мигалку" на имеющих отношение к ней узлах

      if (cols_lo(me) <= 51 .and. cols_hi(me) >= 49) m (50, 50) = 1
      if (cols_lo(me) <= 52 .and. cols_hi(me) >= 50) m (50, 51) = 1
      if (cols_lo(me) <= 53 .and. cols_hi(me) >= 51) m (50, 52) = 1

    end subroutine init_matrix

    ! gfortran считает эту подпрограмму impure из-за доступа в чужую память
    ! ifort и ifx считают её 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

        do j = cols_lo (me), cols_hi (me)
        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

      ! обмениваемся гало с соседями снизу и сверху
      if (me > 1) then
        m2 (:, cols_lo (me)) [me-1] = m2 (:, cols_lo (me))
      end if
      if (me < num_images()) then
        m2 (:, cols_hi (me)) [me+1] = m2 (:, cols_hi (me))
      end if

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

      if (me == num_images()) m2 (:, 0) [1] = m2 (:, cols)
      if (me == 1) m2 (:, cols+1) [num_images()] = m2 (:, 1)
      m2 (0,cols_lo(me):cols_hi(me))       = m2 (rows, cols_lo(me):cols_hi(me))   
      m2 (rows+1, cols_lo(me):cols_hi(me)) = m2 (1, cols_lo(me):cols_hi(me))

    end subroutine process_step

    subroutine output_matrix (m)

      integer (kind=matrix_kind), intent (in) :: m (0:,0:)
      integer :: rows, cols
      integer :: i, j, n
      integer :: outfile

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

      open (file = 'life.txt', newunit=outfile)
      do i = 1, rows
        write (outfile, '(*(A1))') (char (ichar (' ') + &
          m(i, j)*(ichar ('*') - ichar (' '))), j=1, cols)
      end do
      close (outfile)

    end subroutine output_matrix

end program life_mpp

Для трансляции и запуска программ в mpp архитектуре зачастую используются специальные скрипты:

$ caf life_mpp.f90 -o life_mpp -O3 -ftree-vectorize -fopt-info-vec -flto -fopenmp

$ mpiexec -n 4 ./life_mpp

 it's me:           1

 it's me:           2

 it's me:           3

 it's me:           4

           2 сек,    480267000 ячеек/с

Как видим, gfortran ускорил нашу программу ещё на 30% по сравнению с лучшим результатом SMP, достигнув на 4 ядрах производительности 3.84 от последовательной версии. Очевидно, это связано с тем, что даже в SMP системе расшивка параллельных процессов по адресам памяти ускоряет работу с оперативной памятью.

У Intel Fortran:

$ ifort life_mpp.f90 -o life_mpp_ifort -Ofast -coarray

$ ifx life_mpp.f90 -o life_mpp_ifx -Ofast -coarray

$ mpiexec ./life_mpp_ifort

it's me:           1

 it's me:           2

 it's me:           3

 it's me:           4

 it's me:           5

 it's me:           6

 it's me:           7

 it's me:           8

          12 сек,    112970000 ячеек/с

$ mpiexec ./life_mpp_ifx

 it's me:           1

 it's me:           5

 it's me:           6

 it's me:           7

 it's me:           8

 it's me:           3

 it's me:           4

 it's me:           2

          19 сек,     73470000 ячеек/с

ifort в полтора раза хуже, чем его же производительность SMP, про ifx нечего и говорить.

Для полноты картины попробуем ещё забавный режим, запустив ifort с режимами одновременно -coarray и -parallel:

$ mpiexec ./life_mpp_ifort

 it's me:           1

 it's me:           2

 it's me:           3

 it's me:           4

 it's me:           5

 it's me:           6

 it's me:           7

 it's me:           8

          15 сек,     91120000 ячеек/с

Результат предсказуем, SMP и MPP параллелизации в данном случае мешают друг другу, и выполение программы замедляется. Если бы мы запускали комассивную программу на матрице узлов MPP, каждый из которых внутри был бы организован как SMP, тогда такой режим компиляции для ifort имел бы смысл.

Надо отметить, что утилита создания MPP среды mpiexec из комплекта ifort (Intel(R) MPI Library, Version 2021.9  Build 20230307) работает в SMP режиме с ошибкой: вместо задания параметром количества узлов, она всегда создаёт узлы по количеству ядер, а в зависимости от значения параметра дублирует соответствующее число раз эти узлы с теми же самыми номерами, что, разумеется, недопустимо.

Вывод

Мы рассмотрели введение в программирование на Фортране для массивно-параллельных архитектур, а GNU Fortran опять показал своё превосходство над Intel Fortran при использовании языковых конструкций параллельного программирования.

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