Лекция Стэнфорда: проектирование параллельных алгоритмов через примитивы последовательностей

В лекции из цикла CS149 Стэнфордского университета рассматривается фундаментальный сдвиг в подходе к параллельному программированию — переход от низкоуровневого управления потоками к декларативному дата-параллельному мышлению. Профессор объясняет, как сложные, нерегулярные алгоритмы и структуры данных можно эффективно раскладывать на базовые высокопараллельные примитивы, такие как map, fold, scan и segmented scan. Этот подход исключает возникновение скрытых зависимостей данных и лежит в основе архитектуры современных тензорных библиотек и крупномасштабных систем распределенных вычислений.

🌐 Переход к дата-параллельному мышлению и вызовы современных архитектур 0:04

В рамках вводной части занятия лектор обозначает логистические изменения в расписании курса: текущая лекция посвящена преимущественно алгоритмическим андеграундам, а на следующей неделе профессор Кунле Олукотун (Kunle Olukotun) вернет курс к аппаратному обеспечению, начав с разбора фреймворка Apache Spark и вопросов когерентности кэша. Основной акцент сегодняшнего материала смещается с ручного распределения вычислений по потокам на использование богатого набора декларативных примитивов над коллекциями данных.

Масштабы параллелизма в современных вычислительных чипах требуют пересмотра старых парадигм программирования. Профессор напоминает, что один современный графический процессор (GPU) содержит около 163 000 одновременно выполняющихся контекстов. Чтобы полностью утилизировать такую вычислительную мощность и эффективно скрывать задержки доступа к памяти, разработчику необходимо оперировать наборами данных, обеспечивающими параллелизм в режиме не менее 200 000 независимых операций.

Традиционный подход к созданию параллельных программ всегда начинается с поиска и анализа зависимостей в коде. Обнаружив участки, свободные от зависимостей, программист явным образом выделяет под них параллельные потоки. Новая концепция предлагает принципиально иной путь: выразить алгоритм через комбинацию готовых стандартных функций, про которые заранее известно, что они имеют максимально оптимизированную и масштабируемую параллельную реализацию. Ярким примером из повседневной практики программирования является библиотека NumPy в Python или специализированные тензорные языки. При сложении двух больших векторов в NumPy программист не думает о потоках или аппаратных барьерах — он просто вызывает оператор сложения, делегируя параллелизацию низкоуровневой библиотеке.

🔢 Концепция последовательностей и базовые примитивы: Map и Fold 4:38

Для формализации дата-параллельного подхода лектор вводит понятие «последовательности» (sequence) — упорядоченной коллекции элементов. Различные языки программирования и фреймворки по-разному воплощают эту абстракцию:

В C++ существует встроенный тип sequence.
В языке Scala используются списки (lists).
В экосистеме Python применяются DataFrames или массивы NumPy.
В фреймворках машинного обучения PyTorch и TensorFlow эта концепция представлена тензорами.

Ключевое отличие последовательности от стандартного массива заключается в строгом ограничении способов доступа к данным. В обычном массиве программист может в любой момент запросить произвольный элемент по индексу (например, A[43]). Лектор подчеркивает, что именно эта свобода произвольного доступа зачастую порождает скрытые, трудноотслеживаемые зависимости между итерациями циклов, что приводит к некорректной работе параллельного кода. Последовательности же разрешают доступ к своим элементам только через строго определенный набор высокоуровневых операций.

Примитив Map (Отображение)

Первым и наиболее очевидным примитивом является функция высшего порядка map. Она принимает в качестве аргументов входную последовательность и функцию $f$, которая затем применяется к каждому элементу индивидуально, формируя выходную последовательность. С точки зрения типов, если функция $f$ преобразует тип $A$ в тип $B$ ($A \to B$), то map принимает последовательность типа $A$ и возвращает последовательность типа $B$.

В синтаксисе C++ это выглядит несколько громоздко через итераторы и стандартную функцию std::transform, тогда как в функциональных языках типа Haskell запись предельно лаконична. Лектор отмечает, что map является идеальным кандидатом для распараллеливания: по определению, каждый вызов функции $f$ полностью изолирован. Разработчик функции $f$ пишет код для одного элемента, не задумываясь о коллекциях вообще. Простейшая параллельная реализация map силами создателя библиотеки сводится к разбиению последовательности $S$ на $P$ равных частей (по числу процессоров или потоков), последовательной обработке каждого блока и последующей конкатенации результатов.

Примитив Fold (Свертка)

Вторым важнейшим примитивом выступает функция fold, переводящая последовательность элементов в одно скалярное значение. Она последовательно применяет бинарный оператор к элементам коллекции. Примером работы fold с оператором сложения является вычисление суммы всех элементов массива.

На вопрос о возможности параллелизации fold в аудитории возникает дискуссия. Лектор подтверждает, что параллельное вычисление суммы возможно через древовидную структуру промежуточных сумм, но делает важное теоретическое уточнение: произвольная функция свертки гарантированно поддается параллелизации только в том случае, если передаваемый бинарный оператор обладает свойством ассоциативности. При этом коммутативность оператора не является обязательной, если библиотека строго контролирует порядок выполнения вычислений.

Помимо этого, знание семантики примитивов позволяет современным JIT-компиляторам (например, в PyTorch) проводить глубокую оптимизацию кода. Если в программе подряд идут операции map и fold, компилятор автоматически производит слияние (fusion) этих шагов. Вместо двух проходов по огромному массиву в памяти создается один цикл, который выполняет трансформацию элемента и сразу же добавляет его к промежуточному результату свертки, экономя пропускную способность кэша.

🔄 Операция Scan: от последовательного алгоритма к эффективному распараллеливанию 20:01

Операция сканирования (scan), также известная как префиксная сумма, преобразует входную последовательность в выходную последовательность той же длины. Каждый $i$-й элемент на выходе представляет собой результат последовательного применения ассоциативного бинарного оператора ко всем элементам от начала массива до индекса $i$ включительно. Таким образом, последний элемент результирующей последовательности равен результату операции fold, но scan сохраняет также и все промежуточные частичные суммы. Различают инклюзивный (включающий текущий элемент) и эксклюзивный (исключающий текущий элемент) scan.

Последовательная реализация сканирования тривиальна: один цикл проходит по массиву, прибавляя текущее значение к накопленной сумме, что требует $O(n)$ операций выполнения работы. Распараллеливание этой операции представляет собой серьезную алгоритмическую задачу.

Эффективный параллельный алгоритм Блеллока (Blelloch Scan)

Если применить наивный параллельный подход на GPU, где количество потоков равно количеству элементов, затраты на вычисления составят $O(n \log n)$ единиц работы. Рост асимптотической сложности с $O(n)$ до $O(n \log n)$ на больших массивах (например, в миллион элементов) критичен и нежелателен для параллельных систем.

Решением является алгоритм, предложенный Гаем Блеллоком (Guy Blelloch), который сохраняет логарифмическое время выполнения (span) $O(\log n)$, но при этом выполняет всего $O(n)$ совокупной работы. Алгоритм состоит из двух фаз:

Phase 1: Upsweep (Восходящая фаза). Строится сбалансированное дерево промежуточных сумм от листьев к корню, аналогично параллельной свертке. На вершине дерева оказывается общая сумма массива.
Phase 2: Downsweep (Нисходящая фаза). Полученные частичные суммы распределяются и перебазируются обратно по дереву вниз к листьям.

За счет древовидной структуры на первом шаге выполняется $n$ операций, на втором — $n/2$, затем $n/4$ и так далее. Сумма этого телескопического ряда сходится к $O(n)$. Профессор отмечает, что алгоритм Блеллока содержит скрытый коэффициент оверхеда, равный примерно 2, поскольку дерево обходится в обоих направлениях. Кроме того, по мере приближения к корню дерева всё больше вычислительных ядер проистекают в простое. Данный алгоритм студентам предстоит реализовать в качестве разминочного задания в рамках лабораторной работы по CUDA.

Аппаратно-зависимые оптимизации Scan

Стратегия распараллеливания радикально меняется в зависимости от целевой архитектуры:

Двухъядерный процессор общего назначения (CPU). Лектор утверждает, что оптимальным инженерным решением будет разделение массива ровно пополам. Каждое ядро последовательно сканирует свою половину на максимальной скорости работы с памятью ($O(n)$), после чего значение из конца первой половины прибавляется ко всем элементам второй половины параллельным потоком. Это обеспечивает идеальный баланс загрузки без накладных расходов на древовидные структуры.
SIMD-архитектура / CUDA Warp. На уровне одного варпа (32 потока в CUDA, выполняющие инструкции синхронно), алгоритм с избыточной работой $O(n \log n)$ парадоксальным образом оказывается значительно эффективнее алгоритма Блеллока с эффективностью $O(n)$. Нативная реализация сканирования внутри варпа занимает всего 5 инструкций сложения (поскольку $\log_2 32 = 5$). Если переписать этот код под эффективный по работе алгоритм Блеллока, потребуется 5 инструкций на подъем и 5 на спуск (всего 10 циклов). Так как SIMD-ядра все равно простаивают при ветвлениях, избыточная работа в 5 циклов выгоднее, чем 10 циклов со сложной логикой условных переходов.
Иерархический Scan больших массивов на GPU. Реальные промышленные библиотеки (например, стандартные библиотеки NVIDIA) комбинируют эти подходы. Сначала множество варпов параллельно сканируют независимые блоки по 32 элемента (5 циклов). Затем блок из 32 получившихся частичных сумм сканируется еще раз (еще 5 циклов). На финальном этапе каждый поток прибавляет базовое смещение к своему элементу (1 цикл). В итоге массив из 1024 элементов сканируется всего за 11 тактов.

В ответ на вопрос из аудитории о динамическом отключении питания (gating) неактивных SIMD-лайнов для экономии энергии, лектор, ссылаясь на авторитетное мнение Билла Дэлли (Bill Dally) из NVIDIA, подтверждает, что подобные технологии энергосбережения внедряются, однако они дают выигрыш по питанию в пределах коэффициента 1.5 и вторичны по сравнению с алгоритмической оптимизацией утилизации варпов.

📊 Сегментированный сканирующий примитив (Segmented Scan) и разреженные матрицы 45:05

В реальных приложениях разработчикам часто приходится иметь дело не с плоскими массивами, а со структурами типа «последовательность последовательностей» (массивы массивов варьирующейся длины). Примерами служат:

Списки смежности в графах (для каждого узла — свой список ребер).
Физические симуляции частиц (для каждой частицы — список соседей в радиусе взаимодействия).
Текстовая аналитика документов (для каждого документа — список содержащихся в нем слов).

В таких задачах скрыты два уровня параллелизма: внешний (по объектам) и внутренний (по свойствам объектов). Если на GPU распараллелить вычисления только по внешнему уровню (например, если в графе всего 10 000 узлов, но у каждого много ребер), мощностей графического процессора не хватит для полной загрузки.

Примитив segmented scan принимает на вход последовательность последовательностей и выполняет операцию сканирования независимо и параллельно внутри каждого сегмента. На практике такая структура упаковывается в два плоских массива: массив самих данных и массив бинарных флагов, где единица указывает на начало нового сегмента. Алгоритм Блеллока модифицируется: флаги начала сегмента протягиваются по дереву сумм, и если на каком-то этапе восходящее или нисходящее ребро пересекает границу сегмента, передача информации блокируется с помощью простых условных операторов if. Модифицированный сегментированный алгоритм также работает за время $O(\log n)$ при объеме вычислений $O(n)$.

Применение: Умножение разреженной матрицы на вектор (SpMV)

Классическим примером эффективности сегментированного сканирования является операция умножения разреженной матрицы, где большинство элементов равны нулю, на плотный вектор. Для компактного хранения таких матриц используется формат CSR (Compressed Sparse Row). Ненулевые значения строк матрицы вытягиваются в одну плоскую последовательность последовательностей, параллельно формируется массив индексов соответствующих столбцов, а также массив указателей на позиции начала каждой строки.

Реализация алгоритма SpMV через дата-параллельные примитивы включает три шага:

Gather (Сбор данных). На основе массива столбцов из плотного вектора $X$ формируется промежуточный массив значений, по длине равный количеству ненулевых элементов матрицы.
Map (Отображение). Выполняется поэлементное параллельное умножение (умножаются ненулевые значения матрицы на собранные элементы вектора $X$).
Segmented Scan (Сегментированное сложение). Проводится сегментированное сканирование со сложением, где границами сегментов служат маркеры начала строк матрицы.

Последний элемент каждого сегмента на выходе будет содержать итоговую сумму строки. Уровень параллелизма такого алгоритма пропорционален общему числу ненулевых элементов матрицы, а не числу строк, что позволяет эффективно обрабатывать сверхкрупные разреженные структуры данных на GPU.

🔄 Перемещение данных: примитивы Gather, Scatter и эмуляция записи 57:01

Для изменения структуры расположения данных в памяти дата-параллельная модель полагается на два ключевых примитива — gather (сбор) и scatter (распределение). Функция gather принимает массив индексов и исходный массив данных, формируя результирующую последовательность путем извлечения элементов из исходного массива по указанным адресам. Функция scatter выполняет противоположное действие: берет плотный массив данных и распределяет его элементы по разреженным адресам целевого массива большой емкости.

Лектор предупреждает, что эти операции сопряжены с высокими аппаратными издержками, поскольку они инициируют непредсказуемый доступ к памяти, зависящий от обрабатываемых данных. Современные процессоры поддерживают аппаратные SIMD-инструкции для сбора (например, команды Gather в наборе AVX2). Регистр инструкций принимает вектор индексов и базовый указатель, после чего параллельно подтягивает данные.

Однако на физическом уровне выполнение такой инструкции чревато тем, что каждый отдельный SIMD-лайн может спровоцировать независимый промах кэша (cache miss) или вызвать индивидуальную ошибку отсутствия страницы в памяти (page fault). Именно поэтому в ISPC или CUDA последовательный доступ к соседним ячейкам памяти (vector load) предпочтительнее произвольного сбора через gather.

Эмуляция Scatter через Sort и Segmented Scan

В условиях отсутствия аппаратной поддержки бесконфликтной параллельной записи (scatter с возможностью коллизий индексов, как при построении гистограмм) архитекторы ПО применяют изящный алгоритмический трюк. Процесс распределения конфликтующих данных по неуникальным индексам можно эмулировать с помощью сортировки и сегментированного сканирования:

Сортировка (Sort). Входной массив данных параллельно сортируется по значениям целевых индексов. Теперь все элементы, предназначенные для записи по одному адресу, лежат в памяти рядом.
Поиск границ (Map). С помощью простого параллельного сопоставления каждого элемента со своим левым соседом выявляются точки смены индексов, и создается бинарный вектор флагов начала сегментов.
Свертка коллизий (Segmented Scan). Внутри каждого сформированного сегмента проводится сегментированное сканирование с целевым оператором (например, суммированием для гистограммы).
Финальная запись. Результаты агрегации бесконфликтно переносятся в итоговые ячейки памяти.

🌌 Практический кейс: Параллельное построение пространственной сетки для симуляции частиц 1:03:50

В качестве сквозного примера, напрямую соотносящегося с домашним заданием студентов, профессор разбирает задачу построения пространственной сетки для астрофизического моделирования галактик или гидродинамических симуляций частиц. Пространство разбивается на фиксированную сетку (например, 16 ячеек в конфигурации 4 на 4). На входе — миллионы частиц со случайными непрерывными координатами $(x, y)$. Задача — построить структуру данных (сетку списков), которая для каждой ячейки возвращает список идентификаторов содержащихся в ней частиц. Данная структура критически важна для n-body симуляций, чтобы вычислять силы взаимодействия между физическими объектами только в пределах соседних ячеек.

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

Попытка завести независимую копию сетки для каждого воркера с последующим слиянием (merge) хорошо работает на 4-8 ядерных CPU, но абсолютно не масштабируется на GPU с сотнями тысяч потоков, так как аллокация и финальное слияние тысяч сеток поглотят все процессорное время. Изменение оси параллелизма (запуск потоков по ячейкам, где каждая ячейка линейно сканирует миллионы частиц) также неэффективно, поскольку заставляет каждый поток выполнять колоссальный объем лишней работы.

Дата-параллельный алгоритм построения сетки

Проблема нерегулярного параллелизма элегантно разрешается в рамках дата-параллельной парадигмы без использования блокировок за четыре шага:

Шаг 1 (Map): Потоки параллельно берут координаты каждой частицы и вычисляют номер ячейки, в которой она находится. Результат записывается в плоский массив grid cell. На этом этапе коммуникация между потоками отсутствует.
Шаг 2 (Parallel Sort): Массив частиц сортируется по вычисленным номерам ячеек с сохранением исходных ID частиц (группировка через сортировку). Теперь все частицы, попавшие в нулевую ячейку, гарантированно лежат в начале массива, за ними идут частицы первой ячейки и так далее.
Шаг 3 (Map): Запускается параллельный проход, где каждый поток сравнивает номер ячейки текущего элемента с левым соседом. Если номера различаются, поток фиксирует индекс начала новой группы и записывает этот адрес в итоговый массив cell starts.
Шаг 4 (Результат): На выходе формируется компактный массив из 16 элементов, хранящий указатели на начало и конец непрерывных блоков в отсортированном массиве частиц.

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

🛠️ Реальный мир: От библиотек для GPU до экосистемы Apache Spark 1:16:11

В заключительной части лекции профессор демонстрирует, как описанные теоретические абстракции маскируют аппаратную сложность в коммерческих программных комплексах реального мира. На уровне программирования графических чипов компания NVIDIA поставляет высокоуровневую библиотеку Thrust. Изучив документацию Thrust API, разработчик обнаружит те самые функции: thrust::map, thrust::sort, thrust::scan. Библиотека позволяет писать лаконичный параллельный код на C++ без необходимости ручного проектирования сложных CUDA-кернелов, распределения разделяемой памяти и выравнивания варпов.

В области распределенных вычислений на больших кластерах по тем же принципам построена знаменитая система Apache Spark. Ключевая абстракция Spark — RDD (Resilient Distributed Dataset) — по своей сути представляет собой распределенную последовательность данных. Платформа намеренно запрещает программисту произвольный доступ к элементам кластера, предоставляя взамен фиксированный набор декларативных дата-параллельных операторов.

Именно это ограничение позволяет экосистеме Apache Spark брать на себя полную автоматизацию параллелизма по узлам кластера, прозрачно обеспечивать отказоустойчивость, организовывать резервное копирование данных и эффективно оптимизировать распределенные аналитические запросы. Профессор прощается с аудиторией, пожелав успехов в завершении текущей лабораторной работы и анонсируя выдачу нового масштабного задания по CUDA в ближайший понедельник.