Построение множества Мандельброта — классический пример чрезвычайно параллельной задачи (embarrassingly parallel problem).
Вначале мы разберем наивную реализацию, поиграемся с интринсиками (intrinsics) и, не теряя переносимости, заставим компилятор генерировать нам SIMD-инструкции. Далее добавим многопоточность и в заключение обесценим все наши старания несколькими строчками на CUDA.
Предыстория
На первом курсе МФТИ, на Факультете радиотехники и компьютерных технологий, наш преподаватель по программированию, Илья Дединский, познакомил нас с векторными инструкциями x86-64, используя в качестве примера множество Мандельброта и его удивительные свойства. Эта тема сразу меня увлекла, и я захотел углубиться в дальнейшие оптимизации, но в течение семестра мне не хватало ни времени, ни знаний. Спустя год я решил восполнить этот пробел.
Почему множество Мандельброта?
Множество Мандельброта — множество точек c на комплексной плоскости, для которых рекуррентное соотношение zₙ₊₁ = zₙ² + c задаёт ограниченную последовательность при z₀ = 0.
То есть мы берем некоторую точку c на плоскости и проверяем для неё ограниченность последовательности zₙ. На первый взгляд, задача кажется вычислительно невозможной: для каждой точки нам пришлось бы проверять бесконечное количество итераций. К счастью, мы можем ввести ряд упрощений:
Можно доказать, что если на некотором шаге zₙ² > 4 (или |zₙ| > 2), то последовательность гарантированно расходится
Мы ограничимся конечным количеством итераций N: если за N шагов |zₙ| так и не вышел за радиус 2, то считаем, что уже никогда не выйдет.
Но главное преимущество этой задачи — её идеальная распараллеливаемость. Каждый пиксель:
Вычисляется независимо от соседей,
Требует одинаковых операций.
Визуализация
Хоть визуализация множества, строго говоря, не необходима, но, "во первых, это красиво", а, кроме того, с ней можно быстро увидеть ошибки в построении.
Для красивой визуализации придумаем цветовую палитру, которая будет отображать, насколько быстро точка «убежала» за пределы радиуса 2. Если же точка не вышла за радиус за N шагов — присвоим ей черный цвет.
Вот, например, функция, использованная для генерации заставки статьи (используется графическая библиотека SFML версии 2.6)
sf::Color get_color(int iteration, int maxIterations) { if (iteration >= maxIterations) { return sf::Color::Black; } const float colorScale = 255.0f / maxIterations; const float iterNormalized = iteration * colorScale; const sf::Uint8 r = (sf::Uint8)(iterNormalized / 2 + 0); const sf::Uint8 g = (sf::Uint8)(iterNormalized * 2 + 2); const sf::Uint8 b = (sf::Uint8)(iterNormalized * 2 + 5); return sf::Color(r, g, b, 255); }
Разные функции задают разные палитры, это не принципиально, но каждая красива по-своему:
Однако перейдем к оптимизациям.
Глава 1. Однопоточная реализация
Наивная реализация
Приведем алгоритм в псевдокоде:
for each (c_x, c_y) on complex plane: z_x, z_y, iteration = 0.0 while (z_x² + z_y² ≤ 4 AND iteration < 256): xtemp = z_x² - z_y² + c_x # Действительная часть z² + c z_y = 2*z_x*z_y + c_y # Мнимая часть z² + c z_x = xtemp iteration++ set_color_by_iteration(c_x, c_y, iteration)
z_x и z_y — действительная и мнимая части комплексного числа z
c_x и c_y — координаты точки c на комплексной плоскости
Действительную и мнимую часть получаем из: z² = (x + iy)² = x² - y² + 2xyi
На языке C код будет чуть сложнее из-за необходимости явного преобразования между координатами экрана и комплексной плоскости:
void mandelbrot_naive() { for (int screen_y = 0; screen_y < WINDOW_HEIGHT; screen_y++) { for (int screen_x = 0; screen_x < WINDOW_WIDTH; screen_x++) { // [0, 1920] x [0, 1080] -> [-1, 1] x [-1, 1] float c_y = -1.0f + screen_y * (2.0f / WINDOW_HEIGHT); float c_x = -1.0f + screen_x * (2.0f / WINDOW_WIDTH); float z_x = 0.0f, z_x2 = 0.0f, z_y = 0.0f, z_y2 = 0.0f; int iterations = 0; while (z_x2 + z_y2 < MAX_RADIUS_2 && iterations < MAX_ITERATION_DEPTH) { z_y = 2 * z_x * z_y + c_y; z_x = z_x2 - z_y2 + c_x; z_x2 = z_x * z_x; z_y2 = z_y * z_y; iterations++; } sf::Color color = get_color(iterations, 256); draw_pixel(screen_x, screen_y, color); // Рисуем пиксель цветом, // зависящим от к-ва итераций } } }
Ключевые отличия от псевдокода
Отдельно сохраняются z_x² и z_y², полученные для вычисления радиуса, что позволяет:
Избежать повторного вычисления в условии цикла,
Убрать временную переменную xtemp.
Пиксели экрана (screen_x, screen_y) отображаются на область [-1, +1] × [-1, +1] комплексной плоскости.
Этот пример кода и все последующие будут обладать некоторыми упрощениями. Они не учитывают соотношение сторон монитора, не поддерживают масштабирование и смещение. Полную версию примеров со всеми фичами можно найти у меня на GitHub.
Для дальнейшего понимания происходящего, важно хорошо понять пример выше.
Первый бенчмарк (измерение скорости работы)
Методика измерения
Для проверки корректности работы алгоритма будем выводить с помощью SFML картинку и смотреть на неё "глазками" (обычно, когда что-то идёт не так, это крайне легко заметить).
Для большей точности измерения перформанса (производительности) будем замерять время работы алгоритма несколько (N) раз, при этом не выводя ничего на экран, т.к. это занимает много времени (мы же хотим измерить время работы исключетельно алгоритма, а не скорость рисования картинки). Я выбирал N так, чтобы работа программы занимала порядка 5 секунд.
Я использовал утилиту hyperfine, которая позволяет автоматизировать процесс измерений:
Пример вычисления FPS
FPS = 200 кадров / 4.272 секунд ≈ 46.8 fps (кадров в секунду)
dFPS = (0.019 / 4.272) * 46.8fps ≈ 0.2 fps
Итого: (46.8 ± 0.2) fps
Результаты наивной реализации
Тестовая конфигурация
CPU | AMD Ryzen 5 5600H @ 3.3 GHz |
Оперативная память | 16 GB DDR4 @ 3200 MHz |
ОС | Arch Linux (версия ядра 6.13.5-arch1-1) |
Компиляторы | Clang 19.1.7, GCC 14.2.1 |
Разрешение | 1920×1080 (Full HD) |
Глубина итераций | 256 |
Компилятор, оптимизирующие флаги | Результат, FPS |
GCC, -O2 и выше | 7.0 ± 0.1 |
Clang, -O1 и выше | 7.0 ± 0.1 |
Векторизация
Что такое SIMD?
SIMD (Single Instruction, Multiple Data) — принцип компьютерных вычислений, позволяющий обеспечить параллелизм на уровне данных.
Современные CPU на архитектуре x86-64 (например, с поддержкой расширений AVX/AVX2) используют 256-битные регистры (ymm0-ymm15), которые могут хранить 8 чисел типа float.
Что за расширения архитектуры?
Расширение архитектуры процессора — это набор дополнительных инструкций и/или регистров. Обычно, они повышают производительность в конкретных задачах. Векторные расширения (изначально MMX), например, были придуманы для ускорения процессов кодирования/декодирования потоковых аудио- и видеоданных.
Таким образом, сложение двух регистров ymm1 и ymm2 выполняется как:
ymm1 = [a₀, a₁, a₂, a₃, a₄, a₅, a₆, a₇]
ymm2 = [b₀, b₁, b₂, b₃, b₄, b₅, b₆, b₇]
Результат: [a₀+b₀, a₁+b₁, ..., a₇+b₇]
Получается, одной инструкцией мы заменили восемь. В нашем случае это означает, что мы сможем обрабатывать по 8 пикселей за раз, что, в теории, может привести к 8-кратному ускорению.
Пример кода, использующий библиотеку SIMD интринсиков
Пример на псевдокоде (понятнее, читать первым)
c_step_x = 2.0 / WIDTH c_step_y = 2.0 / HEIGHT c_y = -1.0 for screen_y in 0 to HEIGHT: // Векторные переменные будем помечать нижним подчеркиванием _c_x = [-1.0 + 0 * c_step_x, -1.0 + 1 * c_step_x, ..., -1.0 + 7 * c_step_x] for screen_x in 0 to WIDTH step 8: _z_x = [0.0, 0.0, .., 0.0] _z_y = [0.0, 0.0, .., 0.0] _z_x2 = [0.0, 0.0, .., 0.0] _z_y2 = [0.0, 0.0, .., 0.0] _z_xy = [0.0, 0.0, .., 0.0] _iterations = [0, 0, .., 0] for iteration in 0 to MAX_ITERATIONS: _radius_sq = _z_x * _z_x + _z_y * _z_y _cmp_mask = [0.0, 0.0, .., 0.0] for i in 0..7: if radius_sq[i] < MAX_RADIUS_SQ: _cmp_mask[i] = -1.0 else: _cmp_mask[i] = 0.0 if _cmp_mask == [0.0, 0.0, .., 0.0]: break for i in 0..7: _z_x[i] = _z_x2 - _z_y2 + _c_x _z_y[i] = 2 * _z_xy + _c_y _iterations -= _cmp_mask _z_x2 = _z_x * _z_x _z_y2 = _z_y * _z_y _z_xy = _z_x * _z_y for i in 0..7: set_color(screen_x + i, screen_y, iterations[i]) _c_x += [step_x, step_x, ..., step_x] _c_y += step_y
Пример на C++ с использованием библиотеки <x86intrin.h>
#include <x86intrin.h> void mandelbrot_vectorized() { const __m256 _01234567 = _mm256_set_ps(7.0f, 6.0f, 5.0f, 4.0f, 3.0f, 2.0f, 1.0f, 0.0f); // Вместо того, чтобы каждую итерацию цикла заного вычислять c_x и c_y, // будем каждую итерацию увеличивать их на c_step_x и c_step_y. const float c_step_x = (2.0f / WINDOW_WIDTH); // Инициализируем вектор 8ю одинаковыми значениями const __m256 _c_step_x = _mm256_set1_ps(c_step_x); const __m256 _8_c_steps_x = _mm256_set1_ps(c_step_x * 8); const float c_step_y = (2.0f / WINDOW_HEIGHT); const __m256 _c_step_y = _mm256_set1_ps(c_step_y); __m256 _max_radius2 = _mm256_set1_ps(MAX_RADIUS_2); float c_y = -1.0f; __m256 _c_y = _mm256_set1_ps(c_y); for (int screenY = 0; screenY < WINDOW_HEIGHT; screenY++, c_y += c_step_y) { float c_x = -1.0f; __m256 _c_x = _mm256_set1_ps(c_x); // Сдвигаем начальные X-координаты для 8 пикселей: // _c_x = [c_x + 0*step, c_x + 1*step, ..., c_x + 7*step] _c_x = _mm256_add_ps(_c_x, _mm256_mul_ps(_c_step_x, _01234567)); // Обрабатываем по 8 пикселей за раз for (int screenX = 0; screenX < WINDOW_WIDTH; screenX += 8) { // Инициализируем нулями __m256 _z_x = _mm256_setzero_ps(); __m256 _z_y = _mm256_setzero_ps(); __m256 _z_x2 = _mm256_setzero_ps(); __m256 _z_y2 = _mm256_setzero_ps(); __m256 _z_xy = _mm256_setzero_ps(); __m256i _iterations = _mm256_setzero_si256(); for (int iteration = 0; iteration < MAX_ITERATION_DEPTH; iteration++) { __m256 _radius2 = _mm256_add_ps(_z_x2, _z_y2); // Сравниваем радиус² с 4.0, получаем битовую маску: // Для всех 8 float'ов: // -1 (0xFFFFFFFF) если условие выполняется, // 0 (0x00000000) если нет. // Условие - вылетел ли z за максимальный радиус __m256 _cmpMask = _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OQ); // Из каждого int'а маски берем старший бит // 1, если условие выполнилось, // 0, если нет. int mask = _mm256_movemask_ps(_cmpMask); // Если ни для одной точки условие не выполняется (все биты 0), // то прерываем цикл if (mask == 0x00) break; _z_x = _mm256_add_ps(_c_x, _mm256_sub_ps(_z_x2, _z_y2)); _z_y = _mm256_add_ps(_c_y, _mm256_mul_ps(_mm256_set1_ps(2.0f), _z_xy)); // Обновляем счетчики итераций: // Вычтем из вектора итераций вектор маски // Тогда если маска равна 0 (мы вышли за максимальный радиус), // то итерации не меняются // Если маска равна -1 (мы не вышли за максимальный радиус), // то итерации увеличиваются на 1 _iterations = _mm256_sub_epi32(_iterations, _mm256_castps_si256(_cmpMask)); _z_x2 = _mm256_mul_ps(_z_x, _z_x); _z_y2 = _mm256_mul_ps(_z_y, _z_y); _z_xy = _mm256_mul_ps(_z_x, _z_y); } // ... Записываем цвета (пропущено для краткости) // Сдвигаем X-координаты для следующей группы из 8 пикселей _c_x = _mm256_add_ps(_c_x, _8_c_steps_x); } // Обновляем Y-координату для следующей строки _c_y = _mm256_add_ps(_c_y, _c_step_y); } }
Ключевые моменты:
Параллелизм 8:1. Все скалярные операции над
floatзаменены на векторные инструкции с использованием__m256— специального типа данных, который хранит 8 чисел с плавающей точкой.Трюк с маской. Вместо проверки
if (x² + y² < 4)для каждой точки, используется векторное сравнение:__m256 _cmpMask = _mm256_cmp_ps(_radius2, _max_radius2, _CMP_LT_OQ); ... int mask = _mm256_movemask_ps(_cmpMask); if (mask == 0x00) break; ... _iterations = _mm256_sub_epi32(_iterations, _mm256_castps_si256(_cmpMask));Для всех 8 точек одновременно генерируется маска сравнения
_cmpMask(см. псевдокод)-1— условие выполняется,0— условие не выполняется.Из вектора итераций мы вычитаем эту маску. Таким образом, итерация точки инкрементируется только при выполнении условия. Это позволяет обойтись без условных операторов внутри цикла.
Используя инструкцию
_mm256_movemask_ps, мы кладем старшие биты 8ми результатов сравнения в 8-битное число. Если оно равно нулю, то значит все точки вылетели за максимальный радиус, можно переходить к следующей восьмёрке.
Недостатки такого подхода:
Сложность разработки. Написания кода напоминает написание ассемблера — требуется знание интринсиков, особенностей регистров. Постоянно приходится лезть в документацию от Intel.
Совместимость. Код не будет работать на CPU без поддержки расширения AVX, а также на других архитектурах. А если процессор поддерживает AVX-512 (с регистрами размера 512), то не будет использован его полный потенциал.
Это недостатки приводят к тому, что под каждую архитектуру и под каждую ширину векторных инструкций придётся долго и нудно переписывать код.
Результаты ручной векторизации
Компилятор, оптимизирующие флаги | Результат, FPS (относительный прирост) |
GCC, O1 | 42.3 ± 0.1 (х6.0) |
GCC, O2 | 46.4 ± 0.1 (x6.6) |
GCC, O3 | 47.0 ± 0.1 (x6.7) |
Clang, O1 | 42.9 ± 0.1 (х6.1) |
Clang, O2 | 43.2 ± 0.1 (х6.2) |
Clang, O3 | 43.6 ± 0.1 (х6.2) |
Итого, мы получили прирост в 6.7 раз. Конечно, хотелось x8, но не стоит забывать, что время вычисления 8ми пикселей определяется временем самого медленного. Мы хорошо приблизились к теоритическому пределу.
Танцы с компилятором. Разговор глухого и немого
В идеальном мире программист пишет читаемый код, а компилятор волшебным образом превращает его в оптимальный машинный код. Увы, реальность сложнее — современные компиляторы, хоть и умны, часто нуждаются в намёках от программистов. Наша задача — зарефакторить код так, чтобы компилятор увидел векторизацию.
Вместо специальных типов, будем использовать тупо массивы. Обязательно выравняем их в памяти, так как векторные регистры могут лежать только по выравненным адресам. Вместо векторного сложения и умножения, разжуем компилятору сложения массивов независимыми циклами. Также уберём ветвление аналогично прошлому примеру. Теперь осталось только надеятся, что компилятор нас понял.
Пример кода c векторизацией массивами:
const int VEC_SIZE = 8; #define FOR_VEC for (int i = 0; i < VEC_SIZE; ++i) #define ALIGN alignas(VEC_SIZE * 4) void mandelbrot_arrayed() { const float c_step_x = (2.0f / WINDOW_WIDTH); const float c_step_y = (2.0f / WINDOW_HEIGHT); const float vec_c_step_x = c_step_x * VEC_SIZE; float c_y = -1.0f; for (int screenY = 0; screenY < WINDOW_HEIGHT; screenY++) { const float c_x = -1.0f; ALIGN float _c_x[VEC_SIZE] = {}; FOR_VEC _c_x[i] = c_x + c_step_x * i; for (int screenX = 0; screenX < WINDOW_WIDTH; screenX += VEC_SIZE) { ALIGN float _z_x [VEC_SIZE] = {}; ALIGN float _z_y [VEC_SIZE] = {}; ALIGN float _z_xy[VEC_SIZE] = {}; ALIGN float _z_x2[VEC_SIZE] = {}; ALIGN float _z_y2[VEC_SIZE] = {}; ALIGN int _iterations[VEC_SIZE] = {}; for (int it = 0; it < MAX_ITERATION_DEPTH; ++it) { ALIGN float _radius2[VEC_SIZE] = {}; FOR_VEC _radius2[i] = _z_x2[i] + _z_y2[i]; ALIGN int _is_active[VEC_SIZE] = {}; FOR_VEC _is_active[i] = _radius2[i] < MAX_RADIUS_2; // Очевидное ветвление, которое оптимизациями убирает компилятор // По факту: if (_is_active == 0x00) break; bool all_zero = true; FOR_VEC { if (is_active[i] != 0) { all_zero = false; break; } } if (all_zero) { break; } FOR_VEC _z_x[i] = _z_x2[i] - _z_y2[i] + _c_x[i]; FOR_VEC _z_y[i] = 2 * _z_xy[i] + c_y; // Инкрементируем количество итераций только активным пикселям FOR_VEC _iterations[i] += _is_active[i]; FOR_VEC _z_x2[i] = _z_x[i] * _z_x[i]; FOR_VEC _z_y2[i] = _z_y[i] * _z_y[i]; FOR_VEC _z_xy[i] = _z_x[i] * _z_y[i]; } // ... Записываем цвета (пропущено для краткости) FOR_VEC _c_x[i] += vec_c_step_x; } c_y += c_step_y; } }
Получился достаточно читаемый для человека код, мы выполнили нашу часть договора. Осталось проверить, понял ли этот код компилятор.
Результаты стараний компилятора:
Компилятор, оптимизирующие флаги | Результат, FPS (относительный прирост) |
GCC, -O3 | 13.0 ± 0.1 (x1.9) |
Clang, -O1 | 5.2 ± 0.1 (x0.7) |
Clang, -O2 | 45.8 ± 0.1 (x6.5) |
Clang, -O3 | 46.8 ± 0.1 (x6.7) |
Как видно, с GCC общего языка найтись не удалось. Зато Clang меня понял, да ещё как, он даже стал работать быстрее, чем с интринсиками, и почти догнал GCC.
Также я запустил тесты на процессоре с AVX-512 (Intel i5-1135G7), поставив значение переменной VEC_SIZE на 16. И вот результаты:
Реализация, компилятор | Результат, FPS |
Наивная, Clang под -O3 | 5.3 ± 0.1 |
Векторизованная, Clang под -O3 | 66.2 ± 0.5 (x12.5) |
Итого, у нас получилось:
Улучшить читаемость. Код читается несильно сложнее, чем в наивной реализации.
Вернуть переносимость. Для популярных архитектур с помощью препроцессора можно выставить нужные значения переменной VEC_SIZE, а для остальных - оставить единицей, тогда реализация выродится в наивную.
Вернуть справедливость. Теперь всё честно: человек занимается написанием читаемого кода, а компилятор — оптимизацией.
Глава 2. Многопоточная реализация
Ожидаемый прирост
На моём процессоре 6 ядер и 12 потоков. Значит ли это, что мы можем получить прирост в 6 раз? или может в 12? Всё не так однозначно.
Современные процессоры поддерживают технологию гиперпоточности (Hyper-Threading у Intel или SMT у AMD), которая позволяет одному физическому ядру выполнять сразу два программных потока.
Что это значит? Основная идея гиперпоточности заключается в более эффективном использовании вычислительных ресурсов: пока один поток ждёт данные из памяти или по какой-то иной причине простаивает, второй поток может использовать свободные вычислительные блоки ядра.
Например,
Пока поток А записывает данные в память (это очень долгая операция), поток Б может использовать свободные ресурсы: умножать, складывать и т.д.
При полной загрузке процессора гиперпоточность может обеспечить дополнительное ускорение порядка 20–30% по сравнению с использованием только физических ядер (x7-8 в нашем случае). Но стоит готовиться к худшему — в нашем случае, когда узким местом являются арифметические операции, прирост может оказаться ниже, чем для задач, где часто возникают ожидания.
Использование OpenMP
OpenMP (Open Multi-Processing) — открытый стандарт для распараллеливания программ. Даёт описание совокупности директив компилятора (#pragma ...), библиотечных процедур и переменных окружения, которые предназначены для программирования многопоточных приложений на многопроцессорных системах с общей памятью.
В нашем случае OpenMP — это, наверное, самый простой способ распараллелить программу. И делается это буквально в одну строчку:
... #pragma omp parallel for for (int screenY = 0; screenY < WINDOW_HEIGHT; screenY++) { float c_y = -1.0f + c_step_y * screenY; ...
Достаточно убедиться, что каждая итерация цикла независима от всех остальных: для этого я внёс вычисление c_y внутрь цикла. Напомню, что ранее c_y инкрементировалась каждую итерацию цикла, а значит каждая следующая итерация зависила от предыдущей.
Теперь, добавляем директиву #pragma omp parallel for перед циклом. Так, мы говорим компилятору «дружище, раздели этот цикл на 12 равных частей и дай каждому потоку по кусочку». Готово, подключаем библиотеку и тестируем:
Получаем 182 FPS или прирост в ... 4 раза?
Но почему? Чем там эти потоки занимались?
Как оказывается, равное распределение количества итераций цикла по потокам не гарантирует равное распределение вычислительной нагрузки. При статическом распределении (schedule(static), используемом по умолчанию) цикл делится на равные блоки между потоками. Например, для 12 потоков каждый получит по WINDOW_HEIGHT / 12 строк пикселей. Если эти строки попадают в центральную область, поток будет перегружен, в то время как остальные уже завершат работу и будут простаивать.
Динамическое распределение нагрузки
Динамическое планирование решает эту проблему. Параметр schedule(dynamic, N) говорит OpenMP разбить цикл на блоки по N итераций (строк пикселей) и динамически распределять их между потоками. Как только поток завершает обработку своего блока, он сразу получает следующий доступный. Это гарантирует, что все потоки будут постоянно заняты (по крайней мере вначале), даже если отдельные блоки требуют разного времени для вычислений.
Размер блока намеренно увеличен для лучшей читаемости и равен 40.
Это делается также в одну строчку:
... #pragma omp parallel for schedule(dynamic, 8) for (int screenY = 0; screenY < WINDOW_HEIGHT; screenY++) { float c_y = -1.0f * invMagnifier + c_step_y * screenY; ...
Guided планирование
Между статическим и динамическим подходами существует ещё один вариант — guided планирование (schedule(guided, N)). В этом режиме OpenMP автоматически регулирует размер блоков, начиная с крупных и постепенно уменьшая их до минимального значения N. По началу потоки получают крупные блоки итераций, что минимизирует накладные расходы на распределение задач. По мере выполнения цикла размер блоков уменьшается, чтобы даже под конец всем досталось по кусочку и никто не простаивал.
Размер минимального блока равен 8.
Результаты многопоточной оптимизации:
Победителем вышло guided-планирование (хотя, строго говоря, победа в пределах погрешности), несильно обогнав динамическое за счет меньших накладных расходов. В таблице приведем лучшие результаты каждого метода:
Планирование, размер чанка | Результат, FPS (относительный прирост) |
Статическое | 186 ± 7 (x4.0) |
Динамическое, 1, 2, 4 | 385 ± 6 (x8.2) |
Guided, 1, 2, 4 | 388 ± 8 (x8.3) |
Мы получили ускорение в 8.3 раза, а значит, гиперпоточность дала нам, как минимум, дополнительные 38% производительности!
Итоги многопоточной реализации
Вначале мы написали наивную реализацию, осознали множество Мандельброта. Используя библиотеку интринсиков, мы ускорили код в 6.7 раз. Затем, хитро зарефакторив код, мы вернули читаемость и переносимость. А используя OpenMP, нам удалось эффективно распараллелить вычисления на 6 ядер и 12 потоков.
Итого, нам удалось ускорить код в 55.4 раз, с 7 до 388 FPS.
Глава 2.5. Параллелим однопоточную реализацию. ILP
Немного теории
До этого момента мы смотрели на процессор как на простого последовательного исполнителя инструкций. На самом деле любой современный процессор устроен гораздо сложнее. Главные свойства, которые нас сейчас интересуют — это суперскалярность (superscalar) и внеочередное исполнение (out-of-order execution, OoO).
Суперскалярный процессор — это процессор, реализующий форму параллелизма на уровне команд (Instruction-Level Parallelism, ILP) внутри одного ядра. В отличие от простых скалярных процессоров, которые выполняют максимум одну команду за такт, суперскалярный способен отправлять на выполнение сразу несколько инструкций одновременно, распределяя их по разным исполнительным блокам.
Внеочередное исполнение (Out-of-Order execution) — исполнение инструкций не в том строгом порядке, в котором они идут в машинном коде, а по мере готовности данных для них. Это делается для того, чтобы вычислительные блоки не простаивали. Результат внеочередного исполнения обязательно должен быть эквивалентен последовательному.
В качестве примера давайте взглянем на характеристики моего процессора на микроархитектуре Zen 3. Внутри одного ядра у него есть следующие исполнительные блоки:
Исполнительный блок | Количество на ядро |
ALU (Arithmetic Logic Units — целочисленная арифметика) | 4 |
AGU (Address Generation Units — работа с адресами) | 3 |
FPU (Floating Point Units — операции с плавающей точкой) | 4 |
Это значит, что физически процессор способен исполнять одновременно до четырех целочисленных операций и четырех операций с плавающей точкой!
Взгляд назад
Глава 2.5 была написана спустя год после публикации изначальной статьи. В прошлой главе я писал:
При полной загрузке процессора гиперпоточность может обеспечить дополнительное ускорение порядка 20–30% по сравнению с использованием только физических ядер (x7-8 в нашем случае). Но стоит готовиться к худшему — в нашем случае, когда узким местом являются арифметические операции, прирост может оказаться ниже, чем для задач, где часто возникают ожидания.
Однако затем, когда получилось аж 38% прироста производительности от гиперпоточности, я почему-то совершенно не удивился. А стоило!
Действительно, производительность в нашей задаче ограничивает арифметика (такие задачи называют compute-bound, в противовес memory-bound, где узким горлышком являются операции с памятью). Такой большой прирост от Hyper-Threading в чисто математической задаче говорит о том, что ресурсы одного физического ядра сильно недоутилизируются одним потоком.
Почему так происходит? Ответ кроется во внутренних зависимостях кода. Если следующая инструкция в потоке ждет результата работы предыдущей, исполнительные блоки (в нашем случае FPU) стоят без дела, ожидая данных.
Возьмем простой пример:
Рассмотрим псевдоассемблер. Допустим, что все данные уже разложены по регистрам, а в нашем процессоре всего 2 FPU.
a = b / c; // занимает 20 тактов x = y / z; // занимает 20 тактов
Здесь операции независимы. Наш умный процессор это, конечно, заметит, и поэтому эти 2 инструкции будут отправлены на 2 разных FPU одновременно. Вычисление займёт 20 тактов. Оба FPU работали всё время, мы получили 100% производительности ядра
А теперь немного изменим код:
a = b / c; x = a / z; // Внимание: теперь x зависит от a
Что произойдет сейчас? Второе деление не может начаться, пока не закончится первое. Приходится ждать: 20 тактов на первую строчку + 20 тактов на вторую = 40 тактов. Получается, всё время работал только один FPU, а мы потеряли половину производительности на ожидание.
И тут на помощь приходит Hyper-Threading со своим вторым логическим потоком на одно ядро. Пока первый поток ждет окончания своих же вычислений, второй поток забирает простаивающие ресурсы ядра и делает свою независимую математику.
Используя эти знания, вернемся к однопоточной реализации и дооптимизируем её уже в последний раз. Чтобы понять, где именно мы теряем такты, рассмотрим основной векторизованный цикл Мандельброта (именно он крутится почти всё время исполнения нашей программы):
// for (int it = 0; it < MAX_ITERATION_DEPTH; ++it) { ... } .LBB1_4: vaddps ymm12, ymm9, ymm10 vcmpltps ymm12, ymm12, ymm4 vtestps ymm12, ymm12 je .LBB1_6 vsubps ymm9, ymm9, ymm10 vaddps ymm13, ymm9, ymm7 vaddps ymm9, ymm11, ymm11 vaddps ymm11, ymm9, ymm6 vpsubd ymm8, ymm8, ymm12 vmulps ymm9, ymm13, ymm13 vmulps ymm10, ymm11, ymm11 vmulps ymm11, ymm13, ymm11 dec r8d jne .LBB1_4
Инструкции, заканчивающиеся на ps (Packed Single-Precision), исполняются на блоках FPU. Значит ли это, что первые три инструкции (vaddps, vcmpltps и vtestps) будут выполнены одновременно, ведь у нас в запасе целых четыре FPU?
К сожалению, нет. Каждая следующая из этих инструкций использует результат предыдущей: vaddps записывает результат в регистр ymm12, vcmpltps ждет эти данные, читает их и снова пишет в ymm12, а vtestps снова ждет. Получается цепочка зависимости по данным (Data Dependency), которая сводит на нет всё наше количество исполнительных блоков.
Именно здесь на сцену выходит принцип Out-of-Order execution. Ясно, что выполнить одновременно этот зависимый фрагмент мы не сможем. Но ведь процессору никто не мешает временно отложить инструкции, данные для которых еще не готовы, и заглянуть немного вперед. После ветвления (je) у нас идет vsubps ymm9, ymm9, ymm10, которая ждет данные только с предыдущей итерации цикла, а также vpsubd ymm8, ymm8, ymm12. Процессор может начать выполнять их параллельно с первыми тремя. Таким образом, OoO помогает выполнять больше независимых инструкций одновременно.
Однако, как мы уже поняли, даже с OoO исполнением у нашего процессора не получается использовать все свои ресурсы. Так поможем же ему! Почему бы нам не просчитывать в одном цикле сразу два (или больше) независимых вектора пикселей? Пока данные для первого вектора еще в процессе вычисления (и исполнительные блоки простаивают в ожидании), процессор может загрузить эти блоки вычислениями для второго вектора. Фактически, мы отнимаем работу у Hyper-Threading и делаем её прямо в коде на одном логическом потоке.
Как заставить компилятор реализовать такую логику? Очень просто. Вспомним, с чего начиналась наша векторизованная на массивах версия:
const int VEC_SIZE = 8;
Мы можем просто начать увеличивать VEC_SIZE, сделав его 16, 24, 32... Физически процессор всё равно будет использовать 256-битные регистры (шире в AVX2 просто нет), но компилятор сгенерирует код, в котором вычисления для нескольких независимых векторов будут перемешаны. Эта оптимизация называется Instruction Scheduling (планирование инструкций). Компилятор переставит инструкции так, чтобы независимые вычисления чередовались, максимально загружая доступные FPU.
Оптимизация, которую мы сделали вручную, увеличив VEC_SIZE, называется Loop Unrolling (размотка цикла). Мы "размотали" внешний цикл и теперь за каждую итерацию вычисляем в 2 раза больше пикселей.
for (int screenX = 0; screenX < WINDOW_WIDTH; screenX += VEC_SIZE) // было for (int screenX = 0; screenX < WINDOW_WIDTH; screenX += VEC_SIZE * 2) // стало
Замечание про OoO
Справедливости ради, компилятору особо и не нужно изворачиваться с перестановками – процессор сам лучше знает, в каком порядке ему исполнять инструкции для максимальной производительности. Тогда возникает справедливый вопрос: если наш OoO процессор такой умный, зачем мы помогали ему с раскруткой цикла? Разве он не может сам переставить инструкции в удобном порядке?
Контекст, в пределах которого железо может менять порядок команд на лету, ограничен размером буфера переупорядочивания (ReOrder Buffer, ROB). В современных архитектурах он вмещает порядка 200-300 инструкций (например, 256 команд у моего Zen 3). Этого окна с запасом хватает, чтобы переставлять инструкции внутри текущей итерации цикла. Но если мы не развернем код (оставим VEC_SIZE = 8), в ROB поместится только 10-20 итераций внутреннего цикла (из 256 в худшем случае):
const size_t MAX_ITERATION_DEPTH = 256; ... for (int it = 0; it < MAX_ITERATION_DEPTH; ++it) { // Внутри цикла 14 ассемблерных инструкций // 14 * 256 = 3584 инструкций в худшем случае (черные пиксели) }
Поэтому планировщик процессора просто физически не сможет заглянуть так далеко вперед, чтобы «увидеть» независимые инструкции соседних пикселей.
Итоги однопоточной реализации v2
А теперь давайте посмотрим на результаты наших экспериментов с Loop Unrolling'ом:
VEC_SIZE (Unroll) | Результат, FPS (относительный прирост) |
1 | 7.0 ± 0.1 (x1.0) |
8 (1) | 46.8 ± 0.1 (x6.7) |
16 (2) | 78.3 ± 0.2 (x11.2) |
24 (3) | 93.2 ± 0.2 (x13.3) |
32 (4) | 75.3 ± 0.2 (x10.8) |
40 (5) | 51.5 ± 0.1 (x7.4) |
Мы видим ускорение вплоть до Unroll = 3, достигая ускорения в 13.3 раза на одном потоке. Мы ускорили однопоточную реализацию еще в 2 раза! Грубо говоря, до этого мы утилизировали только лишь половину FPU из доступных.
Но почему при unroll 4 и 5 производительность падает? Ответ прост — у нас тупо заканчиваются регистры. На каждую независимую группу из 8 пикселей внутри цикла нам нужно хранить текущие состояния и счетчик итераций, плюс общие константы. При Unroll = 3 нам еще хватает регистров, но при Unroll = 4 — уже нет. Процессору физически не хватает места, и компилятору приходится сбрасывать промежуточные значения на стек. Такое явление называется Register Spilling. Из-за постоянных обращений к памяти мы и видим резкое падение FPS.
И что теперь?
Может показаться, что теперь нам придется переделывать все измерения главы 2, так как мы изменили однопоточную реализацию. К счастью, нет! (хотя я все равно перемерил, на всякий)
В новом однопоточном коде мы достигаем максимума утилизации FPU за счет грамотного использования ILP. В многопоточном же коде мы упирались в тот же самый максимум за счёт Hyper-Threading. Поскольку исполнительные блоки FPU уже работают на пределе в обоих случаях, получить дополнительное ускорение поверх этого физически невозможно
Глава 3. Вычисление на видеокарте. CUDA.
Почему CUDA?
CUDA (Compute Unified Device Architecture) — программно-аппаратная архитектура параллельных вычислений, которая позволяет существенно увеличить вычислительную производительность благодаря использованию графических процессоров фирмы Nvidia.
В частности, мы будем пользоваться компилятором NVCC, который позволяет превращать C/C++ код с расширениями CUDA в машинный код, исполняемый в основном на видеокарте. Одно из преимуществ такого подхода в том, что функции, написанные для вычисления на GPU почти не отличаются от обычных.
Есть множество способов перенести вычисления на видеокарту, но в этой статье мы ограничимся самым простым. Все тесты будут проводиться на видеокарте RTX 2060 super.
Реализация на CUDA
__global__ void mandelbrot_kernel(sf::Uint8* pixels) { // Вычисляем координаты пикселя через идентификаторы блоков и потоков int screen_x = blockIdx.x * 16 + threadIdx.x; int screen_y = blockIdx.y * 16 + threadIdx.y; if (screen_x >= WINDOW_WIDTH || screen_y >= WINDOW_HEIGHT) return; // Далее аналогично наивной реализации: float с_x = -1.0f + screen_x * (2.0f / WINDOW_WIDTH); float c_y = -1.0f + screen_y * (2.0f / WINDOW_HEIGHT); int iterations = 0; float z_x = 0.0f, z_y = 0.0f, z_x2 = 0.0f, z_y2 = 0.0f; while (z_x2 + z_y2 < MAX_RADIUS_2 && iterations < MAX_ITERATION_DEPTH) { z_y = 2 * z_x * z_y + c_y; z_x = z_x2 - z_y2 + c_x; z_x2 = z_x * z_x; z_y2 = z_y * z_y; iterations++; } // ... Записываем цвет пикселя (пропущено) } void mandelbrot_cuda(sf::Uint8* pixels) { size_t size = WINDOW_WIDTH * WINDOW_HEIGHT * 4 * sizeof(sf::Uint8); // Выделяем память на GPU // (Ради краткости примера, память выделяется каждый раз) sf::Uint8* d_pixels = nullptr; cudaMalloc(&d_pixels, size); // Задаем размеры блока и сетки dim3 blockSize(16, 16); dim3 gridSize((WINDOW_WIDTH + 15) / 16, (WINDOW_HEIGHT + 15) / 16); // Запускаем кернел с заданными размерами блока и сетки mandelbrot_kernel<<<gridSize, blockSize>>>(d_pixels); // Копируем данные с GPU cudaMemcpy(pixels, d_pixels, size, cudaMemcpyDeviceToHost); // Освобождаем память cudaFree(d_pixels); }
Ключевые моменты:
Кернел (kernel) — функция, исполняемая на GPU. Чтобы объявить кернел, перед прототипом пишется
__global__.Поток (thread) — базовая единица выполнения. Каждый поток обрабатывает один пиксель, то есть исполняет кернел один раз.
Блок (thread block) — группа потоков. Может иметь от одного до трех измерений. В нашем случае он (16, 16) — двумерный. (Почему именно 16?, P.S. в комментариях к статье более подробный ответ на этот вопрос)
Сетка (grid) — группа блоков. Также может иметь от одного до трех измерений. В нашем случае он выбирается так, чтобы полностью покрыть экран блоками.
Взято из CUDA C++ Programming Guide.
Стрелкой обозначается каждый потокЧтобы эффективно покрыть экран используем трюк с округлением вверх:
dim3 gridSize((WINDOW_WIDTH + 15) / 16, (WINDOW_HEIGHT + 15) / 16);
Формула
screen_x = blockIdx.x * 16 + threadIdx.xвычисляет глобальную координату пикселя для каждого потока, зная номер блока в сеткеblockIdx.xи номер потока в блокеthreadIdx.x.Если ширина или высота экрана в пикселях ровно не делится на блоки по 16, то некоторые потоки в последних блоках будут выходить за границы экрана. Чтобы этого избежать, нужно всегда проверять корректность глобальных координат:
if (screen_x >= WINDOW_WIDTH || screen_y >= WINDOW_HEIGHT) return;Синтаксис вызова кернела немного отличается от вызова функции: помимо агрументов, мы также передаём размер блока и сетки.
mandelbrot_kernel<<<gridSize, blockSize>>>(d_pixels);
Бенчмарк реализации на CUDA
Запускаем тесты — в этот раз процессорный кулер уже не сходит с ума.
Получаем: 950 ± 5 fps (x2.4)?
Маловато. Как оказалось, мы не используем все ресурсы видеокарты. Всё дело в cudaMemcpy — после её вызова GPU простаивает, ожидая завершения передачи. Этого можно избежать, используя асинхронную передачу данных. Тогда мы приблизимся к теоретическому пределу PCI-E 3.0 16x — шины, по которой передаются данные с моей GPU. Это примерно 16 Гбайт/с / (1920x1080x4 байт) ≈ 2000 fps
Но статья не про это. Тем более, для отрисовки на монитор мы зачем-то передаём данные с GPU на CPU и обратно. Для следующего теста будем измерять только время вычислений.
Результат: 6100 ± 50 FPS
Это в 15.7 раз быстрее мультипоточной реализации и в 870 раз быстрее наивной!
Заключение
Мы прошли длинный путь оптимизации вычисления множества Мандельброта, ускорив его вычисление на почти 3 порядка!
Реализация | Результат, FPS (относительный прирост) |
Наивная | 7.0 ± 0.1 |
SIMD - ручная | 47.0 ± 0.1 (x6.7) |
SIMD - переносимая на массивах | 46.8 ± 0.1 (x6.7) |
Многопоточная - статическая | 186 ± 7 (x4.0) |
Многопоточная - динамическая | 385 ± 6 (x8.2) |
Многопоточная - guided | 388 ± 8 (x8.3) |
CUDA | 6100 ± 50 (x15.7) |
Наш случай был, конечно, простейшим, но на этом примере мы рассмотрели путь оптимизаций, который можно проделать для большинства легко-распараллеливаемых задач.
