Современный Фортран представляет собой специализированный язык программирования, предназначенный в основном для написания вычислительных программ для векторно-конвейерных и параллельных архитектур. Эволюция стандартов языка Фортран была рассмотрена в предыдущих статьях – здесь и здесь.
На данный момент действующим стандартом языка Фортран является стандарт 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 обладает более мощным оптимизатором.
Постановка задачи
Напишем простейшую программу для реализации классического клеточного автомата игры "Жизнь". Не будем сейчас париться с вводом и выводом, исходную конфигурацию зададим в самой программе, а результирующую конфигурацию после заданного числа шагов выведем в файл. Нас будут интересовать сами вычислительные шаги клеточного автомата. Эта задача хороша для нас тем, что она позволяет небольшими усилиями достичь любого наперёд заданного объёма чистых (pure) вычислений с массивами произвольных размеров, не вырождаясь в заведомо излишний код, который оптимизатор мог бы выкинуть, обманув наши метрики производительности.
Для тестов, чтобы далеко не ходить, используется компьютер Mac mini с процессором Intel Core i3 @ 3.6 GHz с 4 физическими ядрами. Компиляторы GNU Fortran 12.2.0 и Intel Classic Fortran 2021.8.0 20221120.
0. Последовательная программа
Для начала напишем программу в чисто последовательном стиле. Напишем всё в одном файле, чтобы оптимизатору было легче работать.
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 (0:,0:)
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 -Ofast
Запустим:
$ ./life_seq_g
11 сек, 125172000 ячеек/с
$ ./life_seq
14 сек, 94120000 ячеек/с
125 лямов в секунду у GNU Fortran против 94 лямов у Intel Fortran.
Попробуем запустить автоматический параллелизатор (спасибо@AlexTmp8за замечание в комментариях):
$ gfortran life_seq.f90 -o life_seq_g -O3 -ftree-vectorize -fopt-info-vec -flto -floop-parallelize-all -fopenmp
$ ifort life_seq.f90 ‑o life_seq ‑Ofast ‑parallel
$ ./life_seq_g
11 сек, 124773000 ячеек/с
$ ./life_seq
4 сек, 340690000 ячеек/с
Intel Fortran очень серьёзно прибавил в производительности, в три с половиной раза. GNU Fortran добавил самую малость. Это единственный из наших тестов, где ifort показал преимущество перед gfortran, причём весьма заметное.
Давайте, может, попробуем 32-разрядные целые вместо байтов (с автопараллелизатором)?
integer, parameter :: matrix_kind = 4$ ./life_seq_g
10 сек, 131818000 ячеек/с
$ ./life_seq
6 сек, 212080000 ячеек/с
Как видим, ничего хорошего нам это не дало.
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_g12 сек, 115730000 ячеек/с
$ ./life_mat
Как видим, результат чуть-чуть хуже чисто последовательного алгоритма. Причём если выключить автопараллелизатор, то Intel Fortran почему-то сильно расстраивается:
$ ./life_mat
25 сек, 55580000 ячеек/с
При этом надо ещё отметить, что 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 -Ofast -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 -Ofast -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 поддерживает мощные методы совместной оптимизации раздельно расположенных в исходных файлах единиц компиляции, поэтому для большой программы сравнительный результат мог бы быть другим.
Получается, что компилятор Intel Classic Fortran (ifort) более эффективен тогда, когда нам нужно оттранслировать на SMP много унаследованного последовательного кода, автоматически его распараллелив. GNU Fortran же позволяет генерировать более эффективный код в абсолютном зачёте, но требует для этого некоторой ручной работы по явному указанию параллелизации.
В следующей статье мы рассматриваем средства поддержки массивно-параллельных архитектур, имеющиеся в современном Фортране.
Комментарии (56)
mobi
03.04.2023 14:55А если сравнить "классический" ifort (сейчас он называется "Intel Fortran Compiler Classic") с новым ifx (теперь именно он называется "Intel Fortran Compiler")?
vadimr Автор
03.04.2023 14:55Может быть, я до такого сравнения и дозрею. Но, судя по тому, что в fix даже чисто синтаксически стандарт языка реализовали только несколько месяцев назад, вряд ли там какая-то особо зрелая технология. При этом под macOS его нет, а под Linux, боюсь, как и все нештатные расширения llvm, он будет вставать с большим скрипом.
agalakhov
03.04.2023 14:55+1Там еще точность вычислений проверить стоит. На фортране обязан работать написанный в лоб алгоритм Кэхэна, который, кстати, не работает в лоб на си при дефолтных флагах компиляции.
mobi
03.04.2023 14:55+2Проверил на ноутбуке i5-5200U@2.20GHz (приведен лучший из трех последовательных запусков):
ifort life_seq: 26 сек, 52470000 ячеек/с ifx life_seq: 25 сек, 54270000 ячеек/с ifort life_mat: 43 сек, 32360000 ячеек/с ifx life_mat: 40 сек, 34450000 ячеек/с ifort life_omp: 13 сек, 100740000 ячеек/с ifx life_omp: 17 сек, 80460000 ячеек/с ifort life_con: 13 сек, 100810000 ячеек/с ifx life_con: 17 сек, 80790000 ячеек/с ifort life_con2: 15 сек, 91720000 ячеек/с ifx life_con2: 51 сек, 27180000 ячеек/сВывод: ifort гораздо лучше справляется с распараллеливанием, но ifx может давать лучший результат для последовательного кода (но скорее всего это в пределах погрешности).
vadimr Автор
03.04.2023 14:55Очень интересные результаты, особенно последний. Можно предположить, что проседание в объединённом параллельном цикле может быть связано с неоптимальным выбором порядка доступа к памяти.
Это у вас Windows или Linux?
mobi
03.04.2023 14:55Windows (кстати, если кто будет повторять: размер стека задается при компиляции ключом /F).
exfizik
03.04.2023 14:55Помнится, когда я последний раз плотно работал с Фортраном (давно уже, лет 15 назад), там зарождалась фича coarray для абстракции MPI/OpenMP. Сейчас уже мне не нужно, и времени вникать в это нет, но было бы интересно почитать, чем дело кончилось.
unreal_undead2
03.04.2023 14:55$ gfortran life_seq.f90 -o life_seq_g -O3 -ftree-vectorize -fopt-info-vec -flto
$ ifort life_seq.f90 -o life_seq -O3
Для честности надо было гну компилятору тоже только -O3 оставить. Или заняться подбором ключей для интеловского - отдельное увлекательное занятие )
vadimr Автор
03.04.2023 14:55Интеловский компилятор включает векторные оптимизации по умолчанию в режиме -O3.
unreal_undead2
03.04.2023 14:55Кроме векторизации там ещё много всего есть (сходу не скажу, например, включает ли -O3 IPO/LTO, можно с -ip/ipo поиграться). И в любом случае в последнее время рекомендуется векторизовать явно через omp simd. Сталкивался с тем, что с новой версией icc автовекторизация слетала, думаю в фортране то же самое. Ну и -march всё таки стоит обоим компиляторам сказать, а то сгенерят какой нибудь древний SSE.
PS В современном gcc "Vectorization is enabled at
-O2which is now equivalent to the original-O2 -ftree-vectorize -fvect-cost-model=very-cheap" https://gcc.gnu.org/gcc-12/changes.htmlvadimr Автор
03.04.2023 14:55В целом, у компилятора Intel очень разумно выбраны умолчания режимов оптимизации кода.
По умолчанию оба компилятора генерируют исполняемый модуль, оптимизированный под текущий процессор (-mtune=native). Добавление -march=native или какого-нибудь -march=skylake не ускоряет программы и даже немножко замедляет их (возможно, это связано с какими-нибудь лишними врапперами с библиотечным кодом в таком случае).
Через omp simd можно векторизовать (особенно если хочется использовать GPU и компилятор это поддерживает), но не это является в данном случае нашей целью, так как мы занимаемся параллелизацией.
unreal_undead2
03.04.2023 14:55-mtune влияет только на скедулинг кода, но не даёт использовать инструкции, которых может не быть на другом железе (скажем, AVX или AVX512).
Через omp simd можно векторизовать (особенно если хочется использовать GPU
OMP Offload явно отдельная тема, там кроме векторизации много нюансов ) Но и на CPU выигрыш от omp simd может быть заметный - хотя соглашусь, что это тоже тема для отдельной статьи.
так как мы занимаемся параллелизацией.
А про выставление affinity в следующей серии расскажете? Оно, конечно, на многосокетной системе важнее, но и на клиентской машинке перекидывание потока с ядра на ядро не так безобидно.
AlexTmp8
03.04.2023 14:55+1Вот за ifort было обидно:
ifort life_seq.f90 /O3 /Qparallel
дает ускорение в 6 раз.vadimr Автор
03.04.2023 14:55Справедливое замечание, хотя на 4-ядерном процессоре, конечно, не в 6 раз. Ключ parallel даёт возможность автоматически частично параллелизовать последовательный код, хотя и не так эффективно, как явное указание параллелизма. Отражу это в статье, спасибо! Сравним автоматические параллелизаторы.
AlexTmp8
03.04.2023 14:55+1Да, это на моем 6-ядерном райзене. Я уточнил цифры:
Ryzen 5 ( 6 ядер) - 5,61
i9 (8 ядер) - 5,91как видно ускорение замедляется и узким местом становиться память.
vadimr Автор
03.04.2023 14:55На старом 8-ядерном Xeon 4C E3-1270v2 из второй части статьи – ускорение вообще всего в 2 раза.
AlexTmp8
03.04.2023 14:55в 2 раза это странно. возможно что-то пошло не так.
добавлю Xeonoв:
E-2224 (4 ядра) 3,39
Gold 6140 (18 ядер) 8 раз.
но лучшее время все равно у i9 - 1.9 сек,
6140-й дремал на 1Гц ))
Hemml
03.04.2023 14:55+1Отличная статья, спасибо!
Я бы предложил разбить сетку на куски, помещающиеся в кэш ядра, потом эти куски обрабатывать параллельно, используя SIMD на нижнем уровне, тогда можно выдоить практически всю производительность из процессоров. Подозреваю, автоматический параллелизатор интеловского фортрана так и делает, по этой причине он и дал выигрыш в одном из тестов.
vadimr Автор
03.04.2023 14:55Вы смотрите прямо в корень. Похожая вещь сделана во второй части статьи, и это действительно даёт прирост производительности.
Hemml
03.04.2023 14:55+1Мы занимаемся ГД/МГД-моделированием на сетках и открыли для себя этот эффект довольно давно -- оказалось, что если создать небольшой массив, скопировать туда данные с сетки, обработать их и скопировать результат в результирующую сетку, это может дать прирост производительности в разы, хотя, казалось бы, лишние копирования. Обычно эта "кэш-сетка" у нас одномерная, так как мы делаем прогонки вдоль осей, для вычисления потоков.
Когда-то я ударился в эксперименты и решил написать генератор кода, который на входе получал бы уравнения в матричном виде и на выходе давал бы распаралеленную программу. В итоге идея не взлетела (по причине моего слабого понимания некоторых вещей на тот момент), но позволила мне дешево поэкспериментировать с разными методами распаралеливания, так как оно все равно делалось автоматически. Так я с удивлением узнал, что spatial decomposition + MPI дает бОльшую производительность, чем чистый OpenMP на одном процессоре, несмотря на обмены MPI, за счет эффекта маленькой сетки. В итоге я пришел к тому, что нужно брать одномерные столбцы/строки из трехмерной сетки, объединять их в пучки подходящего размера и отправлять на счет, часть в CPU, часть в GPU, чтобы загрузить вообще всё. Но я не смог тогда справиться с синхронизацией и в результате половина мощностей всегда простаивала. А в те редкие моменты, когда она не простаивала, узел кластера зависал из-за перегрева или из-за какой-то другой причины, связанной с перегрузом :)
vadimr Автор
03.04.2023 14:55Идеально синхронизировать вообще невероятно сложно.
Мы занимаемся ГД/МГД-моделированием на сетках
Это что-то вроде уравнения Навье-Стокса?
Hemml
03.04.2023 14:55Ну примерно, но в астрофизических применениях. Вязкость у нас своеобразная, не такая как в уравнениях Навье-Стокса.
agalakhov
Репутация Intel Fortran как "быстрого компилятора" создавалась еще тогда, когда gfortran как такового не было, был g77 в составе GCC 2.95, который оптимизировал действительно очень так себе. Я в те времена экспериментировал с будущим gfortran в альфа-версии (кажется, он тогда назывался gfc) и сборками ATLAS, и Intel тогда действительно заметно опережал g77 на настройках по умолчанию, но выкручиванием оптимизации можно было к нему приблизиться. Насколько я знаю, у нас я был единственным, кто вообще заморачивался этим вопросом, все остальные просто слепо верили, что Intel быстрее, но не могли сказать, насколько. Подозреваю, с тех времен так же слепо и верят.
thevlad
Еще те кто пишут на фортране свято верят, что по скорости они уделывают C/С++, как бог черепаху, и что вся настоящая быстрая математика написана на фортране. Хотя к примеру в тех же численных библиотеках питона, а конкретнее реализации BLAS(OpenBLAS/MKL/ATLAS), его уже давно закопали.
zzzzzzerg
Если посмотреть OpenBLAS на гитхабе - то там около 30% кода на фортране, а С около 50 - пока не сильно похоже что его давно закопали (https://github.com/xianyi/OpenBLAS)
mobi
Насколько я знаю, это в основе scipy лежат биндинги к фортрановскому коду, и никто там пока ничего менять не планирует.
agalakhov
Такое тоже, хотя тут важнее другой фактор. Фортран - единственный язык, для которого стандарт оговаривает нюансы поведения floating-point. Компилятор C обходится с оптимизацией намного более вольно. Например, C как правило считает, что x + a - x == a, а фортран учитывает пограничный случай с ошибкой округления. Или, например, gcc самовольно заменяет float на double, а фортран, если ему велели работать в одинарной точности, в одинарной и посчитает. В некоторых численных алгоритмах это важно, алгоритм Кэхэна - наглядный пример.
Чтобы отказаться от фортрана наконец, надо в каком-то из современных языков в стандарте специфицировать вычисление выражений с плавающей запятой до последней значащей цифры вне зависимости от оптимизации. Тогда не останется ни одной причины использовать фортран, кроме legacy.
thevlad
Большинство все таки пишет численно стабильные алгоритмы. А так насколько я понимаю в C достаточно скобок набросать, чтобы был нужный порядок.
vadimr Автор
Скобки в Си не задают физический порядок вычислений. Они просто определяют логическую интерпретацию выражения.
thevlad
Нет, скобки в Си однозначно задают "дерево" и порядок вычислений в нем. https://stackoverflow.com/a/17403097
https://stackoverflow.com/a/52788550
vadimr Автор
По вашей же ссылке написано, что порядок вычислений не определён. Это дерево можно обходить в любом порядке. Когда вы пишете, например, a*b, то никакими скобками вы не добьётесь гарантии, чтобы a вычислялось вперёд b.
thevlad
Какая разница в каком порядке обходить дерево, если для подвыражений выполняется отношение "получен результат до". Или вы думаете что обход дерева: ((a1*b1) + (a2*b2)) + (a3*b3))
в порядке:
r1=a1*b1, r2=a2*b2, r3=a3*b3, r4=r1+r2, r5=r4+r3, output r5
и
r3=a3*b3, r2=a2*b2, r1=a1*b1, r4=r1+r2, r5=r4+r3, output r5
будет иметь разный результат? Процессоры так не работают, арифметические операции по определению дают одинаковый результат, для одинаковых входных данных(и иногда плюс минус состояния).
vadimr Автор
Во-первых, у вычислений может быть побочный эффект различной степени глубины или асинхронное действие (то, что в Фортране задаётся описателями impure, pure, simple).
Во-вторых, даже расстановка скобок a*(b*c) в Си не гарантирует, что b*c будет выполнено вперёд a*b.
thevlad
Не очень понимаю, какие могут быть побочные эффекты у арифметических операций в том же Си, чтобы это на что-то влияло.
Расстановка скобок однозначно задает и гарантирует, об этом было собственно написано в этой ссылке https://stackoverflow.com/a/52788550 с приведением пункта стандарта. И это было проверено мною ни однократно в "живой природе", ни одного отклонения от этого поведения зафиксировано не было.
PS: если не верите, можете сходить на какой ни буть godbolt, врубить максимальную оптимизацию -O3, и увидеть как меняется порядок вычислений в зависимости от скобок.
vadimr Автор
Да что вы такое говорите?
Компилятору плевать на скобки на любом уровне оптимизации.
thevlad
О господи, ну нельзя так.
вот выхлоп:
Объясняю, что там происходит:
1) вызов a()
2) вызов b()
3) вызов c()
4) умножение результата b() на результат c()
5) домножение на результат a()
Про неупорядоченность в Си "дерева" вычислений со скобками, это просто какая-то мифология из мира фортран, возможно тридцати летней выдержки, когда стандарта еще особо и не было.
vadimr Автор
Ну а я о чём говорю? Что скобки не задают порядок вычислений.
Если в каком-нибудь Лиспе напишете
(* a (* b c))
то гарантируется, что сначала вычислится a, потом b, потом c, потом правое умножение, потом левое умножение, потому что аргументы функции вычисляются слева направо. А в Си такого нет.
Вот конкретно про С++ расписано.
thevlad
Скобки задают отношение "вычислено до" на "дереве" вычисления, о чем я вам практически в самом начале и написал. И в вашем примере происходит вычисление result_a*(result_b*result_c), именно в соответствии с этим отношением, как и задает выражение скобок. Где здесь противоречие?
А порядок обхода этого дерева, если расставить все скобки не имеет значения, так как это данность современных процессорных архитектур, где арифметические операции условно "атомарны" и зависят лишь от входных аргументов. (что там было во времена мэйнфреймов, на которых гоняли фортран и адаптировали его синтаксис под их железо, мне не ведомо)
vadimr Автор
В моём примере result_a вычисляется до result_b и result_c. Вот то, о чём я пишу. И это важно, когда они являются неатомарными значениями и могут влиять друг на друга. Тем более когда мы обсуждаем параллельное программирование.
А про арифметические операции писал@agalakhov. Я не настолько силён в вещественной арифметике, чтобы в этом вопросе иметь какую-то позицию.
thevlad
Сначала вы утверждали, что в Си нельзя добиться порядка арифметических операций в a*(b*c), и что скобки это логическая конструкция.
Не пишите конструкции в которых больше одного вызова функции в утверждении, если важен порядок их вызова.
vadimr Автор
Да функции здесь вообще не причём.
Предположим, что у нас в какой-нибудь железке переменная timer содержит адрес таймера, считающего такты. Тогда оператор языка Си:
t = *timer - *timer;
присвоит переменной t положительное значение, отрицательное значение или 0?
Ну, допустим, не 0, если timer будет описано, как volatile. Но со знаком (который зависит от порядка вычислений) нет никакой возможности определиться.
Тут для понятности механизма используется вычитание таймера самого из себя, а на практике очень часто при программировании железа встречаются последовательные манипуляции с регистрами, которые важно выполнять в определённом порядке.
Это очень странное свойство для языка низкоуровневого программирования, в качестве которого используется Си. Скорее, этого можно было бы ожидать именно от Фортрана, как языка сильно оптимизируемых компилятором вычислений. Оно не объяснимо ничем, кроме пофигизма K&R, которые вовремя не подумали прописать это в стандарте.
thevlad
timer не может измениться. Если это memory mapped i/o, в каком ни буть микроконтроллере(или ядре и драйвере), или "параллельный" код, то это просто совершенно по другому пишется и реализуется. Фортрана на микроконтроллерах, я к примеру, вообще ни разу не видел, там сплошной Си. (хотя возможно он там и есть)
vadimr Автор
Если вы работаете в модели исполнения кода, когда контекст не может изменяться независимо от хода вычисления значения выражения, то вас вообще не должны волновать вопросы последовательности вычислений.
Напомню, что вы комментируете статью как раз о параллельном коде.
Кстати, с этим связан ещё вопрос взаимовлияния правой и левой частей оператора присваивания, с которым в Си тоже не всё хорошо (например, когда мы присваиваем перекрывающиеся части массива).
thevlad
Вы просто не очень понимаете, о чем идет речь. На Си и плюсах написаны мегатонны параллельного кода. Только он пишется ни так как вам кажется с вашим исключительным опытом в Фортране.
vadimr Автор
Почему вы думаете, что у меня исключительный опыт в Фортране? Я владею несколькими десятками языков программирования и имею 25-летний опыт программирования на Си. У меня есть патент РФ на конструкцию транспьютера, программируемого на Си. Я разрабатывал компилятор диалекта Си. И я уверен, что параллельному программированию на Си я тоже мог бы поучить многих.
Когда вам кажется, что вы открываете собеседникам какие-то неизвестные им глубины понимания, то это не всегда так.
thevlad
Тогда бы вы не написали предыдущий пост про *timer (1) - *timer (2), предполагая что выборка из памяти 1 произойдет раньше чем 2, если вы напишите это на Фортране. Даже если написать на ассемблере с явным указанием порядка операций на современных out-of-order процессорах это не гарантируется. Нужны барьеры или другие средства сериализации.
Вообще "параллельное программирование" это один из самых сложных вопросов в современной разработке и computer science, поэтому я даже начинать не хотел.
vadimr Автор
Вас не затруднило бы процитировать, где я такое предполагал?
Выше я приводил в этом контексте пример на Лиспе, и в Лиспе порядок вычислений однозначно определён.
Никакое внеочередное исполнение (out of order) в процессоре не вызовет перестановку выборок из памяти одним ядром по одному и тому же адресу (в нашем случае *timer) в обратном порядке, просто потому что у одинаковых операций нет причин переупорядочиться в очереди. Ну и естественно, что на отображаемые в память регистры в микроконтроллерах вообще не распространяется оптимизация доступа к памяти.
Ну и не надо вообще смешивать в одну кучу семантику операторов языка программирования и процессорную архитектуру, это запутывает и без того непростой вопрос.
thevlad
Лисп подвержен точно той же проблеме не упорядоченности выборок из памяти, что и любой другой язык.
https://en.wikipedia.org/wiki/Memory_ordering
"On most modern uniprocessors memory operations are not executed in the order specified by the program code."
vadimr Автор
Выборки из памяти в машинном языке могут быть неупорядоченными, а дело транслятора – обеспечить конвенции, принятые в языке высокого уровня, относительно объектов этого языка.
Если речь о Фортране, то задача его компилятора – автоматически расставить всякие там барьеры или точки синхронизации, где они нужны. А если речь о Лиспе, то программы на нём не работают непосредственно с оперативной памятью, они обращаются к атомам, и эти обращения должным образом упорядочиваются средой выполнения. Семантика интерпретации S-выражений в Лиспе гарантирует, что их элементы вычисляются слева-направо, и никак иначе. Вычисление второго операнда просто не может начаться, пока не закончилось вычисление первого (за исключением специальных параллельных форм).
thevlad
Нет никаких конвенций, вы хотя бы ссылки какие ни буть стандарты или еще что-то привили бы. Что на лиспе или чем то еще можно писать многопоточку, и просто делая обычные обращения к памяти получать их линейную сериализацию.(то есть в том порядке в каком они идут в программе или машинном коде)
Это так не работает, барьеры если и могут расставляться автоматически, то расставляются только в том случаи, если объявлена семантика обращения к переменной из разных потоков, или что-то в этом духе. Просто потому что бездумное расставление барьеров будет ронять производительность на порядок.
agalakhov
В физике много алгоритмов с довольно плохой численной стабильностью. Плохо обусловленных задач, для которых пока не разработаны методы регуляризации. В той же квантовой химии интегралы от быстро осциллирующих функций, которые только считать в лоб, потому что сама функция происходит из численного решения столь же гадкой задачи, аналитически тут никак не подступиться. Мы однажды неправильный знак в ответе получить ухитрились чисто из-за набегания погрешности.
И нет, скобок в Си формально недостаточно. Даже разбиения выражения на несколько строк с промежуточными переменными недостаточно. Оптимизатор компилятора имеет полное право код вида
превратить в
double y = a;, и это не будет считаться багом. Несколько версий назад gcc именно так и поступал, сейчас "поумнел", но формально по стандарту он не обязан считать это правильно.thevlad
У меня есть знакомые, которые профессионально пишут вычислительные алгоритмы на си/плюсах. Да иногда порядок важен и другого выхода нет. Обычно это решается, скобочками или разбиением выражения на утверждения(разделенные ";") однозначно задающими порядок.
Сейчас, многие большие вычислительные пакеты, типа ANSYS, адаптируют под GPU, так там порядок вычисления принципиально до конца плохо определен, и с этим мало что можно сделать.
PS: обычно это настраивается флажками оптимизации, отключающими агрессивную оптимизацию на плавучке. Да, мог выкинуть вычисления, не имеющие с точки зрения компилятора побочного эффекта.
agalakhov
Ну тут не столько порядок вычисления, сколько допустимость тех или иных алгебраических преобразований. Например, замена
a-b = -(b-a)все еще валидна, как и заменаa - a = 0, но вот заменаa + b - a - b = 0уже недопустима. Сейчас с ростом популярности NumPy и нейросетей авторы компиляторов Си стали вычислять эти выражения "как в фортране", и это хорошо. Но это в компиляторе, а в стандарте так и стоит "undefined", и, если по вине этого "undefined" спутник на орбиту не выйдет, то виноваты будут точно не авторы компилятора.