Участник:Bannikov Nikolay/Параллельный алгоритм вычисления всех частных сумм
Свойства и структура алгоритма
Общее описание алгоритма
Данный алгоритм предназначен для вычисления всех частных сумм заданного массива данных. Для входного массива:
требуется найти массив:
такой, что каждый его элемент [math]\displaystyle{ S[i] }[/math] равен сумме всех элементов [math]\displaystyle{ X }[/math] от первого до [math]\displaystyle{ i }[/math]-го включительно.
Алгоритм относится к классу алгоритмов сдваивания. Его ключевая идея заключается в том, что на каждом шаге расстояние, на которое передаются и суммируются уже накопленные частичные результаты, удваиваются.
Входные данные представляют собой массив числовых элементов, для которых определена операция сложения, однако алгоритм может быть обобщен на любую ассоциативную бинарную операцию (умножение, нахождение минимума/максимума, логическое И/ИЛИ и.т.д.).
Алгоритм широко используется в параллельных вычислениях.
Математическое описание алгоритма
Пусть дан массив чисел:
Необходимо найти массив
где [math]\displaystyle{ s_i = x_1 + x_2 + ... + x_i }[/math] для всех [math]\displaystyle{ i }[/math] от 1 до n.
Алгоритм выполняется за [math]\displaystyle{ k = log_{2} n }[/math] шагов:
- Каждый процессор [math]\displaystyle{ i }[/math] хранит значение [math]\displaystyle{ s_i = x_i }[/math]
- Для каждого шага [math]\displaystyle{ t }[/math] от 1 до k:
- Вычисляем расстояние [math]\displaystyle{ d = 2^{(t-1)} }[/math]
- Для каждого индекса [math]\displaystyle{ i \gt d }[/math] процессор [math]\displaystyle{ i }[/math] выполняет: [math]\displaystyle{ s_i = s_i + s_{i-d} }[/math]
- Для каждого индекса [math]\displaystyle{ i ≤ d }[/math] значение [math]\displaystyle{ s_i }[/math] остается неизменным
- После [math]\displaystyle{ k }[/math] шагов массив [math]\displaystyle{ S }[/math] содержит все частные суммы
Простой пример для [math]\displaystyle{ n = 4: }[/math]
- Начало: [math]\displaystyle{ [x_1,x_2,x_3,x_4] }[/math]
- Шаг 1 (d=1): [math]\displaystyle{ [x_1, x_1 + x_2, x_2 + x_3, x_3 + x_4] }[/math]
- Шаг 2 (d=2): [math]\displaystyle{ [x_1, x_1 + x_2, x_1 + x_2 + x_3, x_1 + x_2 + x_3 + x_4] }[/math]
Вычислительное ядро алгоритма
Вычислительное ядро параллельного алгоритма вычисления всех частных сумм представляет собой операцию:
где:
- [math]\displaystyle{ s_i }[/math] - текущее значение в процессоре [math]\displaystyle{ i }[/math]
- [math]\displaystyle{ s_{i-d} }[/math] - значение из процессора, находящегося на расстоянии [math]\displaystyle{ d }[/math] слева
- [math]\displaystyle{ d }[/math] - расстояние между процессорами, которое удваивается на каждом шаге
Характеристики вычислительного ядра:
- Простота: Каждая операция - одно сложение
- Регулярность: На каждом шаге выполняется одинаковый тип операций, только меняется расстояние
- Параллельность: Все операции на одном шаге могут выполняться одновременно разными процессорами
- Коммуникационная интенсивность: Основная сложность - в организации быстрого обмена данными между процессорами, а не в самих вычислениях
На первом шаге выполняется [math]\displaystyle{ n - 1 }[/math] сложений, на втором [math]\displaystyle{ n - 2 }[/math] сложений и так далее. Несмотря на то, что общее количество операций [math]\displaystyle{ (n · log_{2} n) }[/math] больше чем в последовательном алгоритме [math]\displaystyle{ n }[/math] операций, время выполнения сокращается с [math]\displaystyle{ n }[/math] шагов до [math]\displaystyle{ log_{2} n }[/math] шагов благодаря параллелизму.
Макроструктура алгоритма
Алгоритм вычисления всех частных сумм не использует в качестве составных частей другие сложные алгоритмы, но может быть представлен как последовательность однотипных макрошагов (или итераций). Каждый макрошаг состоит из двух базовых операций: коммуникационной (передача данных) и вычислительной (суммирование).
Алгоритм представляет собой последовательность [math]\displaystyle{ log_{2} n }[/math] итераций.
Каждая итерация для всех активных процессоров представима следующими двумя этапами:
- Передача данных на расстояние [math]\displaystyle{ d }[/math]
- Сложение: [math]\displaystyle{ s[i] = s[i] + }[/math] полученное_значение
Схема реализации последовательного алгоритма
Для заданного массива чисел [math]\displaystyle{ X = [x_1, x_2,...,x_n] }[/math] требуется вычислить массив частичных сумм [math]\displaystyle{ S = [s_1,s_2,...,s_n] }[/math], где
В последовательном алгоритме частичные суммы можно вычислить весьма просто из рекуррентного соотношения:
Основной вычислительный фрагмент кода на С++:
// Вычисление частичных сумм
s[0] = a[0];
for (int i = 1; i < n; i++)
{
s[i] = s[i - 1] + a[i];
}
- Последовательный алгоритм выполняет [math]\displaystyle{ O(n) }[/math] операций за [math]\displaystyle{ O(n) }[/math] времени на 1 процессоре
- Параллельный алгоритм выполняет [math]\displaystyle{ O(n \cdot \log n) }[/math] операций за [math]\displaystyle{ O(\log n) }[/math] времени на [math]\displaystyle{ n }[/math] процессорах
- При [math]\displaystyle{ n = 1 }[/math] или при выполнении на одном ядре параллельный алгоритм вырождается в последовательный
- Соотношение между способом хранения данных, арифметическими операциями и временем изображено на диаграмме маршрутизации
- Последовательность вычислений перемещается снизу вверх.
- Операции, которые могут выполняться параллельно, показаны на одном и том же вертикальном уровне.
Последовательная сложность алгоритма
Последовательная сложность алгоритма вычисления всех частичных сумм линейна по размеру входных данных: требуется [math]\displaystyle{ n - 1 }[/math] операции сложения, где [math]\displaystyle{ n }[/math] - размер массива и [math]\displaystyle{ O(n) }[/math] обращений к памяти.
Информационный граф
На рис.2 изображён граф алгоритма. В данном случае выполнено суммирование 8 элементов массива.
- Видно, что параллельный алгоритм вычисления всех частных сумм требует [math]\displaystyle{ log_{2} n }[/math] сложений со степенью параллелизма [math]\displaystyle{ n }[/math], если для обеспечения однородности вычислительной схемы на ПЭ (процессорном элементе), для которых нет данных слева, отправлять нулевые значения.
- Каждая операция алгоритма имеет максимально возможную степень параллелизма [math]\displaystyle{ n }[/math]
- Общее число скалярных арифметических операций увеличилось с [math]\displaystyle{ n }[/math] до [math]\displaystyle{ (n · log_{2} n) }[/math], а время выполнения задачи уменьшилось
Таким образом, увеличение степени параллелизма на каждом этапе от 1 до n происходит за счет существенного увеличения числа скалярных арифметических операций.
Ресурс параллелизма алгоритма
Для вычисления всех частных сумм массива порядка [math]\displaystyle{ n }[/math] с помощью алгоритма в его параллельном варианте требуется выполнить [math]\displaystyle{ (\log_{2} n) }[/math] ярусов с убывающим параллелизмом: от [math]\displaystyle{ n - 1 }[/math] до [math]\displaystyle{ 1 }[/math] операции суммирования на ярус. В предположении доступности неограниченного числа необходимых процессоров требуется [math]\displaystyle{ O(n \log n) }[/math] операций при [math]\displaystyle{ p = n }[/math] за [math]\displaystyle{ O(\log n) }[/math] времени. При классификации по высоте ЯПФ алгоритм имеет логарифмическую высоту [math]\displaystyle{ O(\log n) }[/math],а при классификации по ширине ЯПФ алгоритм имеет линейную ширину ярусов [math]\displaystyle{ O(n) }[/math], реализуя массовый координатный параллелизм по индексу процессорных элементов. Максимальное ускорение ограничено [math]\displaystyle{ O(\frac{n}{\log n}) }[/math] из закона Амдала. Наиболее эффективен при [math]\displaystyle{ n ≫ p ≥ \log n }[/math]. Поскольку эффективность алгоритма уменьшается при росте [math]\displaystyle{ n }[/math] и составляет [math]\displaystyle{ O(\frac{1}{\log n}) }[/math].
Входные и выходные данные алгоритма
Входные данные
- Объект: одномерный массив чисел [math]\displaystyle{ X = [x_1,x_2,...,x_n] }[/math]
- Тип: Любой числовой (целые, вещественные, комплексные)
- Размер: [math]\displaystyle{ n }[/math] элементов ([math]\displaystyle{ n ≥ 0 }[/math])
- Объем памяти: [math]\displaystyle{ n × sizeof }[/math](тип)
Выходные данные
- Объект: Массив частичных сумм [math]\displaystyle{ S = [s_1,s_2,...,s_n] }[/math]
- Формула: [math]\displaystyle{ s_i = x_1 + x_2 + ... + x_i }[/math]
- Размер: [math]\displaystyle{ n }[/math] элементов
- Объем памяти: [math]\displaystyle{ n × sizeof }[/math](тип)
Ключевые особенности
- Вход и выход одинакового размера и типа
- Минимум [math]\displaystyle{ 2n }[/math] элементов памяти (вход + выход)
- Можно перезаписать входной массив с целью экономии памяти
- Последний элемент [math]\displaystyle{ s_n }[/math] содержит сумму всех входных чисел
Ограничения
- Переполнение: Риск для целых чисел при больших [math]\displaystyle{ n }[/math]
- Точность: Накопление ошибок для чисел с плавающей точкой
- Память: Максимальный [math]\displaystyle{ n }[/math] ограничен доступной памятью
Свойства алгоритма
Соотношение последовательной и параллельной сложности в предположении доступности неограниченного числа необходимых процессоров, является [math]\displaystyle{ \frac{n}{\log_{2} n} }[/math] отношение линейной к логарифмической, ограничено законом Амдала.
Алгоритм дает существенное теоретическое ускорение, почти линейное с эффективностью [math]\displaystyle{ \frac{1}{\log_{2} n} }[/math].
При этом вычислительная мощность алгоритма, как отношение числа арифметических операций к суммарному объему входных и выходных данных растет как [math]\displaystyle{ \log n }[/math]: при больших [math]\displaystyle{ n }[/math] на один элемент данных приходится все больше операций. Это положительный момент для реализации на ускорителях и сопроцессорах: накладные расходы на пересылку данных относительно стоимости вычислений снижаются.
Для целых чисел алгоритм абсолютно устойчив, получаем точный результат при отсутствии переполнения. Для чисел с плавающей точкой условно устойчив, поскольку накапливается ошибка округления, где погрешность [math]\displaystyle{ O(n · ε) }[/math], где [math]\displaystyle{ ε }[/math] - машинная точность.
Алгоритм сбалансирован, поскольку лишь одна операция сложения, нет умножений/делений и.т.д. Алгоритм полностью детерминирован, поскольку число этапов фиксировано и равно [math]\displaystyle{ \log_{2} n }[/math], набор операций, выполняемых каждым ПЭ на каждом этапе заранее известен и не зависит от входных данных, отсутствуют ветвления по данным, циклы с выходом по точности, использование случайных чисел и.т.п.
В алгоритме нет «массовых рассылок» скалярных величин с большой степенью исхода одной вершины. Это снижает требования к механизмам широковещательной передачи и упрощает эффективный доступ к данным: каждый результат либо остаётся локальным, либо передаётся ограниченному числу соседей по коммуникационной схеме.
В линейной системе координат присутствуют длинные дуги, что требует аккуратного отображения на реальную топологию, но при подходящем выборе топологии эти дуги становятся локальными.
Граф алгоритма регулярно устроен и допускает компактные укладки и конвейеризацию при реализации на специализированных архитектурах.
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Стандартная последовательная реализация использует простую рекуррентную формулу: каждая следующая частичная сумма вычисляется как сумма предыдущей частичной суммы и текущего элемента массива. Это естественное решение с линейной сложностью [math]\displaystyle{ O(n) }[/math], где выполняется ровно [math]\displaystyle{ n - 1 }[/math] операций сложения. Реализация сталкивается с фундаментальным выбором: сохранять исходные данные, создавая отдельный массив [math]\displaystyle{ S }[/math] размера [math]\displaystyle{ n }[/math], тем самым увеличивая использование памяти в 2 раза, либо перезаписывать исходный массив, что позволяет экономить память, однако безвозвратно уничтожает входные данные. Такой выбор критичен при работе с большими массивами, где объем памяти становится ограничивающим фактором.
Для массивов, превышающих размер кэш-памяти процессора, применяется блочная стратегия. Массив разбивается на блоки размером M, где M выбирается таким образом, чтобы блок полностью помещался в кэш L1 или L2. Типичные значения M составляют 16-64 тысяч элементов для 32-битных чисел. Внутри каждого блока последовательно считаются частичные суммы, для блоков хранятся их полные суммы. При переходе к следующему блоку к его элементам добавляется накопленная сумма всех предыдущих блоков. Это улучшает работу с кэш-памятью и даёт удобную основу для параллельных версий.
Говоря про разрядность, часто нет необходимости выполнять все операции в double. На практике разумен вариант: элементы хранятся в float, а накопитель суммы — в double. Это уменьшает ошибки округления при небольших накладных расходах. Для целых данных следует выбирать достаточно широкий тип (обычно 64-битный), чтобы избежать переполнения.
Ресурс параллелизма непосредственно в линейном цикле ограничен (каждая итерация зависит от предыдущей), но за счёт ассоциативности сложения алгоритм можно эквивалентно преобразовать в блочную или другие формы. Такие формы удобны для распараллеливания и реализации на многоядерных системах и кластерах. Векторизация ограничена зависимостью между итерациями, но эффективность улучшается за счёт выравнивания данных, линейного доступа и использования более узких типов (float).
С точки зрения различных платформ:
- на многоядерных SMP-узлах блоки естественно распределяются по потокам
- на векторных и GPU-архитектурах используются деревообразные формы
- на кластерах частичные суммы считаются локально с последующей глобальной коррекцией
Возможные способы и особенности параллельной реализации алгоритма
Ввиду слабой распространённости задачи нахождения частных сумм в библиотеках данный алгоритм обычно не содержится.
Для алгоритма вычисления всех частичных сумм по одномерному массиву естественно использовать иерархический ресурс параллелизма.
Иерархия параллелизма:
- Векторный уровень (SIMD)
Внутри одного потока данные обрабатываются блоками, для которых сначала вычисляются локальные суммы, затем выполняется локальный префикс внутри блока. Это даёт массовый параллелизм на уровне регистров.
- Уровень нитей / потоков (OpenMP)
Массив делится на [math]\displaystyle{ P }[/math] непрерывных сегментов. Каждый поток считает локальные частичные суммы внутри своего сегмента, затем формируется массив из сумм сегментов длины [math]\displaystyle{ P }[/math]. По этому массиву выполняется префиксная сумма (последовательно или параллельно). Полученные "смещения" добавляются ко всем элементам соответствующих сегментов.
Количество операций: примерно [math]\displaystyle{ O(N) }[/math] арифметики + [math]\displaystyle{ O(P) }[/math] для обработки сумм сегментов на потоках.
Объем дополнительных данных - [math]\displaystyle{ O(P) }[/math].
- Уровень процессов (MPI, кластеры)
Аналогично нитям ,каждый MPI-процесс хранит свой участок массива (распределение по непрерывным блокам), локально считает частичные суммы. Далее выполняется коллективная операция MPI_Scan или эквивалентная схема с MPI_Allgather/MPI_Bcast для префиксной суммы по суммам процессов. Процессы добавляют своё смещение ко всем локальным элементам. Объем коммуникаций: по суммам процессов - [math]\displaystyle{ O(P) }[/math] чисел (обычно одно-два double/int64 на процесс), число шагов ~ [math]\displaystyle{ O(\log P) }[/math].
Распределение операций и данных:
Данные удобно распределять по непрерывным блокам (по индексам), что минимизирует пересылки, обеспечивает линейный доступ к памяти, упрощает векторизацию. Операции в основном локальны: суммирование внутри сегмента. Глобальная часть — это только обработка сумм сегментов/процессов (малый объём).
Оценка сложности:
Для массива длины [math]\displaystyle{ N }[/math], при [math]\displaystyle{ P }[/math] параллельных единицах (нити/процессы):
- Локальные операции: ≈ [math]\displaystyle{ \frac{N}{P} }[/math] сложений на каждую нить/процесс
- Глобальные операции по суммам блоков: ≈ [math]\displaystyle{ P }[/math] сложений в сумме, либо [math]\displaystyle{ \log P }[/math] шагов в дереве
- Объем пересылок (MPI): один-два числа на процесс, не зависящий от [math]\displaystyle{ N }[/math].
Таким образом, вычислительная сложность остаётся [math]\displaystyle{ O(N) }[/math], а масштабируемость в основном ограничивается коммуникацией по [math]\displaystyle{ P }[/math] и пропускной способностью памяти.
Использование конкретных технологий:
Параллельная реализация алгоритма частичных сумм строится как двух–трёхступенчатое сканирование: локальное (внутри блока/процесса), глобальное по блокам/процессам и последующая корректировка, что одинаково хорошо ложится на модели OpenMP, MPI, CUDA и схемы векторизации.
Результаты прогонов
Алгоритм реализовывался с помощью технологии параллельного программирования MPI, запускался на суперкомпьютере "Ломоносов-2".
| Количество MPI-процессов(p) | T | S(p) | E(p) |
|---|---|---|---|
| 1 | 0,00001 | 1 | 1 |
| 2 | 0,0000253 | 0,395 | 0,198 |
| 4 | 0,0000273 | 0,366 | 0,092 |
| 8 | 0,0000396 | 0,250 | 0,032 |
| 16 | 0,0001 | 0,010 | 0,00063 |
| 32 | 0,000224 | 0,045 | 0,0014 |
| 64 | 0,000216 | 0,046 | 0,00072 |
| 128 | 0,000227 | 0,044 | 0,000344 |
| 256 | 0,000102 | 0,098 | 0,000382 |
Выводы для классов архитектур
В целом, для алгоритма частичных сумм на всех архитектурах ключевой фактор — организация памяти (линейный доступ, блочная/локальная работа, минимизация пересылок). Конкретные технологии (OpenMP, MPI, CUDA, директивы SIMD) в основном «оборачивают» одну и ту же двух–трёхэтапную схему: локальное сканирование, префикс по агрегированным суммам и корректировка результатов.
Литература
- А.В.Старченко, В.Н.Берцун. Методы параллельных вычислений. Изд-во Том. ун-та, 2013. - 223 с. (Серия "Суперкомпьютерное образование")