Финальные проекты: битонная сортировка и Мандельброт

1. Сортировка без ветвлений: битонная сеть

Быстрая сортировка на GPU буксует: её сравнения зависят от данных — куда пойдёт элемент, известно только в момент сравнения, а ветвящиеся треды внутри варпа выполняются по очереди (SIMT, C02). Хочется алгоритм, где расписание сравнений фиксировано заранее и не зависит от данных вовсе. Такие алгоритмы называются сортирующими сетями, и самая знаменитая из них — битонная (Кен Бэтчер, 1968).

Битонной называется последовательность, которая сначала не убывает, а затем не возрастает (или наоборот): например, {3, 5, 8, 9, 7, 4, 2, 1}. Ключевое свойство: если к битонной последовательности длины n применить «полуочистку» — сравнить элемент i с элементом i + n/2 и поставить меньший влево, — обе половины станут битонными, причём все элементы левой не превосходят всех элементов правой. Рекурсивно повторяя, за log₂ n раундов получаем отсортированную последовательность. А битонную последовательность из произвольной строят так же: сортируют половины во встречных направлениях. Итого сеть состоит из log₂ n «этапов», этап k — из log₂ k проходов, и на каждом проходе все n/2 сравнений независимы — по сравнению на тред:

Фиолетовые сравнения ставят меньший элемент вверх, оранжевые — вниз: встречные направления первых этапов и создают битонность. Для n = 8 сеть насчитывает 6 проходов и 24 сравнения; для миллиона элементов — 210 проходов по полмиллиона одновременных сравнений каждый. Больше сравнений, чем у quicksort (O(n log² n) против O(n log n)), — но все они ложатся на треды без единого зависимого ветвления. Классический размен этой траектории: больше работы, зато вся параллельная.

2. Ядро: вся магия в i ^ j

Каждый проход сети — отдельный запуск ядра: параметры k (размер сливаемых кусков) и j (текущий шаг сравнения) приходят с хоста, тред i решает, с кем сравниваться, одной битовой операцией:

__global__ void bitonic_step(float* a, int j, int k)
{
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int ixj = i ^ j;          // партнёр: тот же индекс с перевёрнутым битом j

    if (ixj > i) {                     // пару обрабатывает только «левый» тред
        bool up = ((i & k) == 0);      // бит k задаёт направление куска
        if (up ? (a[i] > a[ixj]) : (a[i] < a[ixj])) {
            float t = a[i]; a[i] = a[ixj]; a[ixj] = t;   // обмен
        }
    }
}

// хост: расписание сети — два вложенных цикла
for (int k = 2; k <= n; k <<= 1)          // этапы: 2, 4, 8, …, n
    for (int j = k >> 1; j > 0; j >>= 1) // проходы этапа: k/2, k/4, …, 1
        bitonic_step<<<n / 1024, 1024>>>(dev_a, j, k);

Разберём три строки, в которых всё дело. i ^ j — исключающее ИЛИ: у числа i переворачивается ровно тот бит, который установлен в j (j — всегда степень двойки). Пары «i ↔ i с перевёрнутым битом» — это и есть сравнения прохода: при j = 1 соседи, при j = 4 — элементы на расстоянии 4. Условие ixj > i оставляет из каждой пары одного работника — иначе оба треда попытались бы обменять одну пару дважды. Наконец, i & k проверяет бит k: у первой половины каждого куска длины k он нулевой (сортируем вверх), у второй — единичный (вниз) — так чередуются направления, создающие битонность.

Почему между проходами не нужен __syncthreads? Потому что синхронизация здесь — сам запуск ядра: следующий bitonic_step не начнётся, пока не завершился предыдущий (очередь устройства последовательна). Это стандартный способ глобальной синхронизации «всех со всеми», которой внутри одного запуска не существует (C05): барьер между блоками — это граница между ядрами.

Для сравнения масштабов: соседний проект из студенческого репозитория — чёт-нечётная сортировка (odd-even transposition) — устроена ещё проще: n проходов «сравни соседей через одного», тоже без зависимых ветвлений. Но проходов у неё n, а не log² n: на миллионе элементов это миллион запусков ядра против двухсот десяти. Битонная сеть — редкий случай, когда более хитрая математика окупается тысячекратно.

3. Мандельброт: по треду на пиксель

Второй проект — из мира красоты. Множество Мандельброта — точки c комплексной плоскости, для которых итерация z ← z² + c, начатая с нуля, не убегает в бесконечность. Рисуют его так: каждому пикселю — своя точка c; крутим итерацию, пока |z| ≤ 2, но не больше MAX_ITER раз; число сделанных итераций превращается в цвет. Пиксели полностью независимы — чистейший map на двумерном гриде из C07:

#define W 900
#define H 600
#define MAX_ITER 200

__global__ void mandelbrot(unsigned char* iters)
{
    int px = blockIdx.x * blockDim.x + threadIdx.x;
    int py = blockIdx.y * blockDim.y + threadIdx.y;
    if (px >= W || py >= H) return;               // защита из C04, обе оси

    // пиксель → точка комплексной плоскости [-2.2, 0.8] × [-1, 1]
    float cr = -2.2f + 3.0f * px / W;
    float ci = -1.0f + 2.0f * py / H;

    float zr = 0, zi = 0;
    int k = 0;
    while (k < MAX_ITER && zr * zr + zi * zi <= 4.0f) {
        float t = zr * zr - zi * zi + cr;         // z = z² + c, вручную
        zi = 2 * zr * zi + ci;                    // по действительной и мнимой части
        zr = t;
        k++;
    }
    iters[py * W + px] = (unsigned char)(k * 255 / MAX_ITER);
}

Хост запускает ядро сеткой dim3 grid((W+15)/16, (H+15)/16), block(16,16), забирает массив и сохраняет картинку. Оригинальный студенческий проект показывал её через OpenCV — тяжёлую библиотеку компьютерного зрения, которую ради одной картинки ставить жалко. Обойдёмся без неё: формат PPM — это картинка в виде «заголовок + голые байты», её пишет любой C++ без зависимостей, а открывает любой просмотрщик (и Paint, и GIMP; конвертация в PNG — один клик):

// хост, после cudaMemcpy результата в unsigned char img[W*H]:
FILE* f = fopen("mandelbrot.ppm", "wb");
fprintf(f, "P6\n%d %d\n255\n", W, H);        // заголовок: формат, размер, максимум
for (int i = 0; i < W * H; i++) {
    unsigned char v = img[i];
    unsigned char rgb[3] = { v, (unsigned char)(v / 2), (unsigned char)(128 - v / 2) };
    fwrite(rgb, 1, 3, f);                   // нехитрая палитра из яркости
}
fclose(f);

Пуант проекта — в цифрах. Пикселей 540 тысяч, у каждого до 200 итераций: на CPU однопоточный цикл делает это за секунды, GPU — за единицы миллисекунд, потому что пиксели считаются не по очереди, а по 540 тысяч сразу (насколько позволяет карта). Причём здесь вычислений на каждый байт результата приходится очень много — сотни операций на один записанный байт: по критериям рентабельности из C06 это идеальная задача для видеокарты. Любопытная деталь: соседние пиксели у границы множества делают сильно разное число итераций — треды варпа расходятся по while, и варп работает со скоростью самого медленного (SIMT, C02). Фрактал — живая иллюстрация дивергенции варпа; сгладить её — первое упражнение для любопытных.

4. Куда дальше: галерея настоящих проектов

Все проекты ниже — из студенческого репозитория CUDA-projects; каждый — кандидат на вашу собственную работу по траектории:

• Трассировка лучей (ray tracing). По треду на пиксель, но вместо итераций — луч, летящий сквозь сцену из сфер: пересечения, тени, отражения. Тот же паттерн, что Мандельброт, только математика богаче.
• Сглаживание SSAA (суперсэмплинг). Кадр рендерится в учетверённом разрешении и усредняется: по треду на итоговый пиксель, внутри — мини-редукция 4×4 из C05.
• Классификация по расстоянию Махаланобиса. Каждому объекту — тред, считающий расстояния до центров классов через ковариационные матрицы: перемножения матриц из C07 в статистических декорациях.
• Ранг матрицы. Метод Гаусса на GPU: параллелится обработка строк внутри шага, но сами шаги последовательны — поучительный пример алгоритма, параллельного лишь наполовину (вспомните классификацию из C04).

Выбирая проект, гоняйте его через три вопроса траектории: что здесь map, где прячется редукция и окупит ли объём вычислений копирование данных. Если на все три есть ответ — CUDA вам по плечу.

5. Тренажёр: сеть на живых числах

Контрольные вопросы

• Почему quicksort плохо ложится на GPU, а битонная сеть — хорошо?
  У quicksort сравнения и перемещения зависят от данных — треды варпа расходятся по веткам и выполняются по очереди, а само разбиение последовательно по своей природе. У битонной сети расписание сравнений фиксировано заранее и не зависит от данных: каждый проход — n/2 независимых сравнений, по одному на тред, без зависимых ветвлений.
• Что делают строки ixj = i ^ j и if (ixj > i) в ядре?
  i ^ j переворачивает в номере треда бит j (j — степень двойки) — так тред находит партнёра по сравнению на расстоянии j. Условие ixj > i оставляет из пары одного работника: без него оба треда пары выполнили бы обмен дважды, вернув всё как было (или испортив данные при гонке).
• Как синхронизируются проходы сети, если __syncthreads между блоками не существует?
  Границей между запусками ядра: каждый проход — отдельный bitonic_step, а очередь устройства последовательна — следующее ядро не начнётся до завершения предыдущего. Запуск ядра — стандартный способ глобальной синхронизации всех тредов со всеми.
• Почему сеть требует длину — степень двойки и что делать с другими длинами?
  Пары строятся битовой операцией i ^ j, и структура сети определена для 2^m элементов; на рваной длине партнёры выпадают за границы массива. Решение — дополнить массив до ближайшей степени двойки значениями «плюс бесконечность» (FLT_MAX): паддинг всплывёт в конец и отрежется при выгрузке результата.
• Почему Мандельброт — идеальная задача для GPU по критериям C06?
  Пиксели независимы (чистый map), а плотность вычислений огромна: до сотен итераций на один записанный байт результата, при этом входных данных нет вообще — координаты вычисляются из номера треда. Перевозка ничтожна, работа массивна — обратная ситуация к антипримеру «+1 всем элементам».
• Откуда у Мандельброта дивергенция варпов и где она сильнее всего?
  Число итераций while зависит от точки: внутри множества — все MAX_ITER, снаружи — единицы. Треды одного варпа с разным числом итераций выполняются со скоростью самого медленного. Сильнее всего это у границы множества, где соседние пиксели требуют кардинально разного числа итераций.

Источники

1. Batcher K. E. Sorting networks and their applications // Proceedings of the AFIPS Spring Joint Computer Conference. — 1968. — P. 307–314.
2. Bitonic sorter // Wikipedia : [сайт]. — URL: https://en.wikipedia.org/wiki/Bitonic_sorter (дата обращения: 09.07.2026).
3. CUDA C++ Programming Guide // NVIDIA Docs : [сайт]. — URL: https://docs.nvidia.com/cuda/cuda-c-programming-guide/ (дата обращения: 09.07.2026).
4. Netpbm format (PPM) // Wikipedia : [сайт]. — URL: https://en.wikipedia.org/wiki/Netpbm (дата обращения: 09.07.2026).