Векторизация и оптимизация

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

Векторизация — преобразование цикла так, чтобы одна SIMD-инструкция процессора обрабатывала сразу несколько элементов массива (например 4 или 8 чисел) за один такт, многократно ускоряя обработку данных.

Что такое SIMD и зачем он

Современные процессоры умеют выполнять одну операцию над несколькими данными разом — это называется SIMD (Single Instruction, Multiple Data). Векторные регистры (SSE, AVX, AVX-512) вмещают 4, 8 или 16 чисел, и одна инструкция складывает или умножает все эти числа параллельно. Для численного кода, где одни и те же операции применяются к миллионам элементов массивов, это прямой путь к ускорению в 4-16 раз. Векторизация — это превращение скалярного цикла (по одному элементу за итерацию) в векторный (по целой пачке). Хорошая новость: современные компиляторы Fortran векторизуют автоматически — но только если код написан так, что компилятор может доказать безопасность преобразования. Задача программиста — не мешать ему и давать подсказки.

Что мешает векторизации

Компилятор векторизует цикл, только если уверен, что итерации независимы и порядок их выполнения не важен. Несколько вещей этому препятствуют. Зависимости по данным: если итерация использует результат предыдущей (a(i) = a(i-1) + 1), векторизовать нельзя — порядок существенен. Возможный алиасинг: если два указателя могут ссылаться на одну память, компилятор не знает, безопасно ли менять порядок. Вызовы непрозрачных процедур внутри цикла: компилятор не знает, что делает процедура, и не рискует. Сложное управление: непредсказуемые ветвления и выходы из цикла мешают. Понимая эти препятствия, код пишут так, чтобы их избегать.

! Векторизуется легко: итерации независимы, простая арифметика
real :: a(n), b(n), c(n)
do i = 1, n
  c(i) = a(i) * b(i) + 1.0
end do

! НЕ векторизуется: зависимость a(i) от a(i-1)
do i = 2, n
  a(i) = a(i-1) + b(i)     ! рекуррентность — порядок важен
end do

do concurrent: явное обещание независимости

Fortran 2008 ввёл конструкцию do concurrent — способ сообщить компилятору, что итерации цикла независимы и могут выполняться в любом порядке (в том числе параллельно/векторно). Это снимает с компилятора бремя доказательства независимости — вы берёте его на себя.

real :: a(n), b(n), c(n)
integer :: i

! программист гарантирует независимость итераций:
do concurrent (i = 1:n)
  c(i) = a(i) * b(i) + sqrt(abs(a(i)))
end do

do concurrent — это контракт: вы обещаете, что итерации не зависят друг от друга, не имеют побочных эффектов, конфликтующих между итерациями, и могут идти в любом порядке. Если обещание ложно (например, скрытая зависимость), поведение не определено. Зато компилятор получает свободу векторизовать и распараллелить агрессивно. Это мост между обычными циклами и явным параллелизмом, к которому мы вернёмся в разделе про параллелизм.

Чистые процедуры помогают оптимизатору

Атрибут pure у процедуры — обещание, что она не имеет побочных эффектов: не меняет глобальное состояние, не делает ввод-вывод, её результат зависит только от аргументов. Это ценнейшая информация для компилятора. pure-функцию можно безопасно вызывать внутри векторизуемого цикла, переупорядочивать вызовы, выносить из цикла инвариантные. Объявляя процедуры pure везде, где это возможно, вы открываете оптимизатору больше возможностей.

pure function activation(x) result(y)
  real, intent(in) :: x
  real :: y
  y = 1.0 / (1.0 + exp(-x))    ! сигмоида, без побочных эффектов
end function activation

! компилятор может векторизовать цикл с вызовом pure-функции:
do concurrent (i = 1:n)
  out(i) = activation(in(i))
end do

Ещё строже — elemental-процедуры: pure-функции, написанные для скаляра, но автоматически применимые к массивам поэлементно. Объявив activation как elemental, можно писать просто out = activation(in) для всего массива, и компилятор сгенерирует векторный код. Это вершина выразительности: математическая запись плюс автоматическая векторизация.

Флаги компилятора и измерение

Оптимизация не включена по умолчанию. Компилятор оптимизирует на запрошенном уровне: -O0 (без оптимизации, для отладки), -O2 (разумный баланс, обычный выбор), -O3 (агрессивно, включает векторизацию), плюс целевые флаги вроде -march=native (использовать все инструкции конкретного процессора, включая AVX). Без этих флагов даже идеально написанный код не векторизуется.

# типичная сборка релиза с gfortran:
gfortran -O3 -march=native -funroll-loops program.f90 -o program

# отчёт о векторизации (что удалось/не удалось):
gfortran -O3 -fopt-info-vec-all program.f90 -o program

Но ключевое правило оптимизации — измерять, а не гадать. Интуиция о том, где «тормозит», почти всегда обманывает. Сначала профилируют (gprof, perf, встроенные таймеры через system_clock), находят реальные горячие точки, и только их оптимизируют. Оптимизировать код, на который приходится 1% времени, — пустая трата усилий. Fortran даёт встроенный таймер высокого разрешения system_clock для измерения участков.

integer(8) :: t1, t2, rate
call system_clock(t1, rate)
! ... измеряемый участок ...
call system_clock(t2)
print *, "Время, с:", real(t2 - t1) / real(rate)

Как работает под капотом

Векторизация — это работа бэкенда компилятора. Он анализирует цикл, строит граф зависимостей и, если итерации независимы, генерирует SIMD-инструкции: вместо скалярного addss (сложить одно число) — векторный addps/vaddps (сложить 4 или 8 чисел). Цикл при этом «разворачивается» по ширине вектора: основной цикл обрабатывает элементы пачками по 8, а оставшийся «хвост» (если n не делится на 8) — скалярно. do concurrent и pure работают как обещания, которые позволяют анализатору пропустить дорогостоящее или невозможное доказательство независимости и сразу применить векторизацию. Флаг -march=native сообщает, какие векторные инструкции доступны на целевой машине: без него компилятор консервативно использует лишь базовый набор (например, SSE2), а с ним — широкие AVX-512, обрабатывающие 8 чисел real(8) разом. Понимание этой механики объясняет, почему одни и те же строки кода работают в разы быстрее при правильных флагах и почему рекуррентные зависимости — непреодолимый барьер для векторизации.

Память, выравнивание и препятствия векторизации

Чтобы помогать компилятору векторизовать, полезно понимать дополнительные факторы за пределами зависимостей по данным. Первый — выравнивание (alignment): SIMD-инструкции эффективнее всего работают с данными, чей адрес кратен ширине вектора (например, 32 или 64 байтам). Невыровненные загрузки медленнее или требуют дополнительных инструкций. Компилятор часто сам заботится о выравнивании выделяемых массивов, но в сложных случаях это может стать узким местом, о котором сообщают отчёты векторизации. Второй фактор — предполагаемый алиасинг: если процедура принимает два указателя или два target-массива, компилятор обязан допустить, что они могут перекрываться, и потому не рискует векторизовать. Объявление аргументов как обычных (не pointer) ассумированных массивов, а ещё лучше с атрибутом contiguous, сообщает компилятору, что память непрерывна и не пересекается, открывая путь к векторизации.

subroutine saxpy(n, a, x, y)
  integer, intent(in) :: n
  real, intent(in) :: a
  real, intent(in),    contiguous :: x(:)   ! непрерывный, не алиасит
  real, intent(inout), contiguous :: y(:)
  integer :: i
  do concurrent (i = 1:n)
    y(i) = a * x(i) + y(i)        ! легко векторизуется
  end do
end subroutine saxpy

Атрибут contiguous — обещание, что массив занимает непрерывный участок памяти (а не является, например, разреженным срезом с шагом). Это и помогает векторизации, и позволяет компилятору генерировать более простой код доступа. Третий нюанс — хвостовые итерации: когда число элементов не кратно ширине вектора, основной цикл идёт пачками, а остаток обрабатывается скалярно. Для очень коротких массивов накладные расходы на «разгон» вектора и хвост могут съесть выигрыш — векторизация по-настоящему окупается на длинных циклах в сотни и тысячи итераций. Понимание этих деталей объясняет, почему один цикл векторизуется отлично, а похожий — нет, и помогает писать код, дружелюбный к SIMD: непрерывные массивы, независимые итерации, достаточная длина, минимум ветвлений внутри.

Дисциплина оптимизации: измеряй, не угадывай

Завершая тему производительности, стоит сформулировать инженерную дисциплину, которая важнее любого отдельного приёма. Оптимизация без измерения — это гадание, и оно почти всегда промахивается. Знаменитое правило гласит: преждевременная оптимизация — корень многих бед. Интуиция о том, где программа проводит время, обманчива: программисты регулярно вылизывают код, на который приходится доля процента времени, и не замечают настоящее узкое место. Правильный цикл оптимизации строг и эмпиричен. Сначала — профилирование: запустить программу под профайлером (gprof, perf, Intel VTune) или расставить таймеры system_clock и узнать фактическое распределение времени. Затем — найти горячие точки, те немногие участки (обычно 10-20% кода), где тратится 80-90% времени. И только их — оптимизировать, причём измеряя эффект каждого изменения: стало ли реально быстрее, а не «должно было».

Этот подход уберегает от двух типичных ошибок. Первая — оптимизировать не то: усложнять и запутывать код ради ускорения участка, который и так незаметен в общем времени. Вторая — оптимизировать вслепую: вносить «улучшение», которое на деле ничего не даёт или даже вредит (современные компиляторы и процессоры сложны, и «очевидные» ускорения порой оказываются замедлениями). Помните и про иерархию выигрышей: смена алгоритма (например, с O(n²) на O(n log n)) обычно даёт несравнимо больше, чем микрооптимизация констант, — поэтому сначала убедитесь, что алгоритм правильный, и лишь потом точите его реализацию. Связка из этого и предыдущих уроков складывается в цельную картину высокопроизводительного Fortran: выбирайте хороший алгоритм, пишите в массивном стиле с правильным порядком обхода памяти, опирайтесь на встроенные функции и BLAS/LAPACK, включайте оптимизацию компилятора и векторизацию, а решения об их применении принимайте на основе профилирования, а не догадок. Именно эта дисциплина, а не знание отдельных трюков, превращает корректный код в по-настоящему быстрый — и отличает профессионального вычислителя от любителя, который оптимизирует наугад.

Частые ошибки

• Ожидать векторизацию без флагов оптимизации. На -O0/-O1 векторизации обычно нет. Для скорости нужен хотя бы -O2/-O3 и часто -march=native.
• Скрытые зависимости в цикле. Рекуррентность (a(i)=a(i-1)+...) блокирует векторизацию принципиально; переформулируйте алгоритм, если это горячий путь.
• Непрозрачные вызовы внутри цикла. Обычная (не pure) процедура мешает векторизации. Помечайте чистые процедуры pure/elemental.
• Оптимизировать без профилирования. Усилия уходят не туда. Сначала измерьте (system_clock, профайлер), найдите горячие точки, потом оптимизируйте.
• Лгать в do concurrent. Если итерации на деле зависимы, do concurrent даёт неопределённое поведение. Обещание независимости должно быть истинным.

Итоги

• Векторизация (SIMD) обрабатывает несколько элементов одной инструкцией — ускорение в 4-16 раз на массивных вычислениях.
• Компилятор векторизует автоматически, если итерации независимы; рекуррентности и непрозрачные вызовы этому мешают.
• do concurrent явно обещает независимость итераций; pure/elemental снимают барьеры для оптимизатора.
• Нужны флаги оптимизации (-O3, -march=native) — без них векторизации нет.
• Оптимизируйте по данным профилирования, а не по интуиции; измеряйте через system_clock и профайлеры.