Разгон Мандельброта: SIMD с бубнами, OpenMP и CUDA

Построение множества Мандельброта — классический пример чрезвычайно параллельной задачи (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, которая позволяет автоматизировать процесс измерений:

Результаты наивной реализации

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

Что такое SIMD?

SIMD (Single Instruction, Multiple Data) — принцип компьютерных вычислений, позволяющий обеспечить параллелизм на уровне данных.

Современные CPU на архитектуре x86-64 (например, с поддержкой расширений AVX/AVX2) используют 256-битные регистры (ymm0-ymm15), которые могут хранить 8 чисел типа float.

Таким образом, сложение двух регистров 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

#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);
    }
}

Ключевые моменты:

1. Параллелизм 8:1. Все скалярные операции над float заменены на векторные инструкции с использованием __m256 — специального типа данных, который хранит 8 чисел с плавающей точкой.

2. Трюк с маской. Вместо проверки 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-битное число. Если оно равно нулю, то значит все точки вылетели за максимальный радиус, можно переходить к следующей восьмёрке.

Недостатки такого подхода

1. Сложность разработки. Написания кода напоминает написание ассемблера — требуется знание интринсиков, особенностей регистров. Постоянно приходится лезть в документацию от Intel.

2. Совместимость. Код не будет работать на CPU без поддержки расширения AVX, а также на других архитектурах. А если процессор поддерживает AVX-512 (с регистрами размера 512), то не будет использован его полный потенциал.

Это недостатки приводят к тому, что под каждую архитектуру и под каждую ширину векторных инструкций придётся долго и нудно переписывать код.

Результаты ручной векторизации

Итого, мы получили прирост в 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;
    }
}

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

Результаты стараний компилятора:

Как видно, с GCC общего языка найтись не удалось. Зато Clang меня понял, да ещё как, он обогнал производительность своей ручной оптимизации и почти догнал GCC.

Также я запустил тесты на процессоре с AVX-512 (Intel i5-1135G7), поставив значение переменной VEC_SIZE на 16. И вот результаты:

Итого, у нас получилось:

1. Улучшить читаемость. Код читается несильно сложнее, чем в наивной реализации.

2. Вернуть переносимость. Для популярных архитектур с помощью препроцессора можно выставить нужные значения переменной VEC_SIZE, а для остальных - оставить единицей, тогда реализация выродится в наивную.

3. Вернуть справедливость. Теперь всё честно: человек занимается написанием читаемого кода, а компилятор — оптимизацией.

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

Глава 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 итераций (строк пикселей) и динамически распределять их между потоками. Как только поток завершает обработку своего блока, он сразу получает следующий доступный. Это гарантирует, что все потоки будут постоянно заняты (по крайней мере вначале), даже если отдельные блоки требуют разного времени для вычислений.

Это делается также в одну строчку:

...    
#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. По началу потоки получают крупные блоки итераций, что минимизирует накладные расходы на распределение задач. По мере выполнения цикла размер блоков уменьшается, чтобы даже под конец всем досталось по кусочку и никто не простаивал.

Результаты многопоточной оптимизации:

Победителем вышло guided-планирование (хотя, строго говоря, победа в пределах погрешности), несильно обогнав динамическое за счет меньших накладных расходов. В таблице приведем лучшие результаты каждого метода:

Мы получили ускорение в 8.3 раза, а значит, гиперпоточность дала нам дополнительные 38% производительности!

Итоги вычисления на процессоре

Вначале мы написали наивную реализацию, осознали множество Мандельброта. Используя библиотеку интринсиков, мы ускорили код в 6.7 раз. Затем, хитро зарефакторив код, мы вернули читаемость и переносимость. А используя OpenMP, нам удалось эффективно распараллелить вычисления на 6 ядер и 12 потоков.

Итого, нам удалось ускорить код в 55.4 раз, с 7 до 388 FPS.

Глава 3. Вычисление на видеокарте. CUDA.

Почему CUDA?

CUDA (Compute Unified Device Architecture) — программно-аппаратная архитектура параллельных вычислений, которая позволяет существенно увеличить вычислительную производительность благодаря использованию графических процессоров фирмы Nvidia.

В частности, мы будем пользоваться компилятором NVCC, который позволяет превращать C/C++ код с расширениями CUDA в машинный код, исполняемый в основном на видеокарте. Одно из преимуществ такого подхода в том, что функции, написанные для вычисления на GPU почти не отличаются от обычных.

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

Реализация на 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?)

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

  

  • Чтобы эффективно покрыть экран используем трюк с округлением вверх:

    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 16x — шины, по которой передаются данные с GPU. Это примерно 16 Гбайт/с / (1920x1080x4 байт) ≈ 2000 fps.

Но статья не про это. Тем более, для отрисовки на монитор мы зачем-то передаём данные с GPU на CPU и обратно. Для следующего теста будем измерять только время вычислений.

Результат: 6100 ± 50 FPS
Это в 15.7 раз быстрее мультипоточной реализации и в 870 раз быстрее наивной!

Заключение

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

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