Префиксная сумма

Языковые и библиотечные реализации сканирования могут быть инклюзивными или эксклюзивными. Инклюзивное сканирование включает в себя ввод x_i при вычислении вывода y_i (${\displaystyle y_{i}=\bigoplus _{j=0}^{i}x_{j}}$), в то время как эксклюзивное сканирование его не включает (${\displaystyle y_{i}=\bigoplus _{j=0}^{i-1}x_{j}}$). В последнем случае реализации либо оставляют y₀ неопределенным, либо принимают особое значение «x_-1», с помощью которого можно начать сканирование. Эксклюзивные сканы являются более общими в том смысле, что инклюзивное сканирование всегда может быть реализовано в терминах эксклюзивного сканирования (путем последующего комбинирования x_i с y_i), но эксклюзивное сканирование не всегда может быть реализовано в терминах инклюзивного сканирования, как в случае максимума префиксной суммы.

В следующей таблице приведены примеры функций инклюзивного и эксклюзивного сканирования, предоставляемых несколькими языками программирования и библиотеками:

Языки/библиотеки	Инклюзивное сканирование	Эксклюзивное сканирование
Haskell	scanl1	scanl
MPI	MPI_Scan	MPI_Exscan
C++	std::inclusive_scan	std::exclusive_scan
Scala	scan
Rust	scan

Схематическое представление вычисления 16-входной параллельной префиксной суммы

Hillis and Steele представляют следующий алгоритм параллельной суммы префиксов:[8]

for ${\displaystyle i\gets 0}$ to ${\displaystyle \lceil \log _{2}n\rceil -1}$ do

for ${\displaystyle j\gets 0}$ to ${\displaystyle n-1}$ do in parallel

if ${\displaystyle j<2^{i}}$ then

${\displaystyle x_{j}^{i+1}\gets x_{j}^{i}}$

else

${\displaystyle x_{j}^{i+1}\gets x_{j}^{i}+x_{j-2^{i}}^{i}}$

Обозначение ${\displaystyle x_{j}^{i}}$ означает значение j-го элемента массива x на шаге i.

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

Схематическое представление эффективного вычисления 16-входной параллельной префиксной суммы

Эффективный параллельный алгоритм вычисления префиксной суммы может быть реализован следующим способом:[3][9][10]

1. Вычислить суммы последовательных пар элементов, в которых первый элемент пары имеет четный индекс: z₀ = x₀ + x₁, z₁ = x₂ + x₃ и т. д.
2. Рекурсивно вычислить префиксную сумму w₀, w₁, w₂, ..., последовательности z₀, z₁, z₂, ...
3. Выразить каждый член последовательности y₀, y₁, y₂, ... как сумму до двух членов этих промежуточных последовательностей: y₀ = x₀, y₁ = z₀, y₂ = z₀ + x₂, y₃ = w₀, и т. д. Каждое значение y_i, кроме первого, вычисляется либо копированием с позиции вдвое меньшей, чем в последовательности w, либо добавлением значения последовательности x к предыдущему значению y-1_i.

Если входная последовательность имеет размер n, то рекурсия продолжается до глубины O(log n), которая также ограничена временем параллельного выполнения этого алгоритма. Количество операций алгоритма равно O(n), и их можно реализовать на абстрактной параллельной ЭВМ с общей памятью (PRAM) с количеством процессоров O(n/log n) без какого-либо асимптотического замедления путем назначения нескольких индексов каждому процессору в вариантах алгоритма, для которых больше элементов, чем процессоров.[3]

Каждый из предыдущих алгоритмов выполняется за O(log n). Тем не менее, первый делает ровно log₂ n шагов, а второй требует 2 log₂ n − 2 шага. Для приведенных примеров с 16 входами Алгоритм 1 является 12-параллельным (49 единиц работы, разделенных на 4), а Алгоритм 2 — только 4-параллельным (26 единиц работы, разделенных на 6). Однако Алгоритм 2 эффективен с точки зрения работы он выполняет только постоянный коэффициент (2) объема работы, требуемого последовательным алгоритмом, а Алгоритм 1 неэффективен он выполняет асимптотически больше работы (логарифмический коэффициент), чем требуется последовательно. Следовательно, Алгоритм 1 предпочтительнее, когда возможна реализация большого количества параллельных процессов, иначе имеет преимущество Алгоритм 2.

Параллельные алгоритмы для префиксных сумм часто могут быть обобщены на другие ассоциативные двоичные операции сканирования,[3][4], и они также могут быть эффективно вычислены на современном параллельном оборудовании, таком как GPU (Графический процессор).[11] Идея создания в аппаратном обеспечении функционального блока, предназначенного для вычисления многопараметрической префиксной суммы, была запатентована Узи Вишкиным.[12]

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

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

Следующий алгоритм предполагает модель машины с разделяемой памятью; все обрабатывающие элементы PE (от англ. processing elements) имеют доступ к одной и той же памяти. Вариант этого алгоритма реализован в многоядерной стандартной библиотеке шаблонов (MCSTL)[13][14], параллельной реализации стандартной библиотеки шаблонов C++, которая предоставляет адаптированные версии для параллельных вычислений различных алгоритмов.

Чтобы одновременно вычислить префиксную сумму из ${\displaystyle n}$ элементов данных с ${\displaystyle p}$ элементами обработки, данные делятся на ${\displaystyle p+1}$ блоков, каждый из которых содержит ${\displaystyle {\frac {n}{p+1}}}$ элементов (для простоты мы будем считать, что ${\displaystyle n}$ делится на ${\displaystyle p+1}$). Обратите внимание, что, хотя алгоритм делит данные на ${\displaystyle p+1}$ блоков, параллельно работают только ${\displaystyle p}$ обрабатывающих элементов.

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

Смещенные ${\displaystyle p}$, которые хранятся в последней позиции каждого блока, накапливаются в собственной префиксной сумме и сохраняются в последующих позициях. Если ${\displaystyle p}$ мало, то последовательный алгоритм выполняется достаточно быстро, для больших ${\displaystyle p}$ этот шаг может выполняться параллельно.

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

function prefix_sum(elements) {
    n := size(elements)
    p := number of processing elements
    prefix_sum := [0...0] of size n

    do parallel i = 0 to p-1 {
        // i := индекс текущего PE
        from j = i * n / (p+1) to (i+1) * n / (p+1) - 1 do {
            // Здесь хранится префиксная сумма локальных блоков
            store_prefix_sum_with_offset_in(elements, 0, prefix_sum)
        }
    }

    x = 0

    for i = 1 to p {
        x +=  prefix_sum[i * n / (p+1) - 1] // Построение префиксной суммы по первым p блокам
        prefix_sum[i * n / (p+1)] = x       // Сохрание результатов для использования в качестве смещений во втором цикле
    }

    do parallel i = 1 to p {
        // i := индекс текущего PE
        from j = i * n / (p+1) to (i+1) * n / (p+1) - 1 do {
            offset := prefix_sum[i * n / (p+1)]
            // Рассчитывание префиксной суммы в качестве смещения суммы предыдущих блоков
            store_prefix_sum_with_offset_in(elements, offset, prefix_sum)
        }
    }

    return prefix_sum
}

Алгоритм префиксной суммы гиперкуба[15] хорошо адаптирован для платформ с распределенной памятью и использует обмен сообщениями между элементами обработки. Предполагается, что в алгоритме участвуют ${\displaystyle p=2^{d}}$ PE , равного числу углов в ${\displaystyle d}$-мерном гиперкубе.

Различные гиперкубы для разного количества узлов

На протяжении всего алгоритма каждый PE рассматривается как угол в гипотетическом гиперкубе со знанием общей префиксной суммы ${\displaystyle \sigma }$, а также префиксной суммы ${\displaystyle x}$ всех элементов до самого себя (в соответствии с упорядоченными индексами среди PE), каждая в свой гиперкуб.

• Алгоритм начинается с предположения, что каждый PE является единственным углом нульмерного гиперкуба и, следовательно, ${\displaystyle \sigma }$ и ${\displaystyle x}$ равны сумме локального префикса его собственных элементов.
• Алгоритм продолжается объединением гиперкубов, смежных по одному измерению. Во время каждого объединения ${\displaystyle \sigma }$ обменивается и суммируется между двумя гиперкубами, что сохраняет инвариант того, что все PE в углах этого нового гиперкуба хранят общую префиксную сумму объединенного гиперкуба в своей переменной ${\displaystyle \sigma }$. Однако только гиперкуб, содержащий PE с более высокими индексами, добавляет это ${\displaystyle \sigma }$ к своей локальной переменной ${\displaystyle x}$, сохраняя инвариант. ${\displaystyle x}$ хранит только значение суммы префикса всех элементов в PE с индексами, меньшими или равными их собственному индексу.

В ${\displaystyle d}$-мерном гиперкубе с ${\displaystyle 2^{d}}$ PE по углам алгоритм должен повторяться ${\displaystyle d}$ раз, чтобы ${\displaystyle 2^{d}}$ нульмерных гиперкуба были объединены в один ${\displaystyle d}$-мерный гиперкуб. Предполагая модель дуплексной связи, в которой ${\displaystyle \sigma }$ двух смежных PE в разных гиперкубах можно обменивать в обоих направлениях за один шаг связи, это означает, что ${\displaystyle d=\log _{2}p}$ запусков связи.

i := Index of own processor element (PE)
m := prefix sum of local elements of this PE
d := number of dimensions of the hyper cube

x = m;     // Инвариант: префиксная сумма PE в текущем вложенном кубе
σ = m;     // Инвариант: префиксная сумма всех элементов в текущем вложенном кубе

for (k=0; k <= d-1; k++){
    y = σ @ PE(i xor 2^k)  // Получение общей суммы префикса противоположного вложенного куба по измерению k
    σ = σ + y              // Суммирование префиксных сумм обоих вложенных кубов

    if (i & 2^k){
        x = x + y  // Суммирование префиксной суммы из другого вложенного куба только в том случае, если этот PE является индексом с более     высоким индексом
    }
}

Конвейерный алгоритм двоичного дерева[16] — это еще один алгоритм для платформ с распределенной памятью, который особенно хорошо подходит для сообщений большого размера.

Как и алгоритм гиперкуба, он предполагает особую коммуникационную структуру. PE гипотетически расположены в двоичном дереве (например, в дереве Фибоначчи) с инфиксной нумерацией в соответствии с их индексом в PE. Связь в таком дереве всегда происходит между родительским и дочерним узлами.

Инфиксная нумерация гарантирует, что для любого данного PE_j индексы всех узлов достижимы его левым поддеревом ${\displaystyle \mathbb {[l...j-1]} }$ меньше ${\displaystyle j}$, а индексы ${\displaystyle \mathbb {[j+1...r]} }$ всех узлов в правом поддереве больше ${\displaystyle j}$. Индекс родителя больше, чем любой из индексов в поддереве PE_j, если PE_j — левый потомок, и меньше, если PE_j — правый потомок. Это допускает следующие рассуждения:

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

• Локальная префиксная сумма ${\displaystyle \mathbb {\oplus [l..j-1]} }$ левого поддерева должна быть объединена для вычисления суммы локального префикса PE_j ${\displaystyle \mathbb {\oplus [l..j]} }$.
• Локальная префиксная сумма ${\displaystyle \mathbb {\oplus [j+1..r]} }$ правого поддерева должна быть объединена для вычисления суммы локального префикса PE_h более высокого уровня, которая достигается на пути, содержащем левое дочернее соединение (это означает, что ${\displaystyle h>j}$).
• Общая сумма префикса ${\displaystyle \mathbb {\oplus [0..j]} }$ для PE_j необходима для вычисления общей суммы префикса в правом поддереве (например, ${\displaystyle \mathbb {\oplus [0..j..r]} }$ для самого высокого индексного узла в поддереве).
• PE_j должен включать общую префиксную сумму ${\displaystyle \mathbb {\oplus [0..l-1]} }$ первого узла более высокого порядка, который достигается по восходящему пути, включающему соединение правых дочерних элементов, чтобы вычислить его общую префиксную сумму.

Обратите внимание на различие между поддерево-локальной и общей префиксной суммой. Точки два, три и четыре могут привести к тому, что они сформируют круговую зависимость, но это не так. PE нижних уровней могут потребовать, чтобы общая префиксная сумма PE верхних уровней вычисляла их общую префиксную сумму, но PE верхних уровней требуют только локальную префиксную сумму поддерева для вычисления их общей префиксной суммы. Корневой узел как узел самого высокого уровня требует только локальную префиксную сумму своего левого поддерева для вычисления своей собственной префиксной суммы. Каждый PE на пути от PE₀ к корневому PE требует только локальную префиксную сумму своего левого поддерева для вычисления своей собственной префиксной суммы, в то время как каждому узлу на пути от PE_p-1 (последнего PE) до PE_root требуется общая префиксная сумма его родителя, чтобы вычислить его собственную общую префиксную сумму.

Это приводит к двухфазному алгоритму:

Восходящая фаза
Распространение локальной префиксной суммы поддерева ${\displaystyle \mathbb {\oplus [l..j..r]} }$ его родителю для каждого PE_j.

Нисходящая фаза
Распространение эксклюзивной (эксклюзивный PE_j, а также PE в его левом поддереве) суммарной префиксной суммы ${\displaystyle \mathbb {\oplus [0..l-1]} }$ всех PE с более низкими индексами, которые не включены в адресуемое поддерево PE_j, на PE более низких уровней левого дочернего поддерева PE_j. Распространение инклюзивной префиксной суммы ⊕ [0…j] на правое дочернее поддерево PE_j.

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

k := number of packets in a message m of a PE
m @ {left, right, parent, this} := // Сообщения в разные PE

x = m @ this

// Восходящая фаза- рассчитывание локальной префиксной суммы поддерева
for j=0 to k-1: // Конвейерная передача: для каждого пакета сообщения
    if hasLeftChild:
        blocking receive m[j] @ left // Замена локального m[j] полученным m[j]
        // Совокупная инклюзивная локальная префиксная сумма из PE с более низкими индексами
        x[j] = m[j] ⨁ x[j] 

    if hasRightChild:
        blocking receive m[j] @ right
        // Мы не объединяем m[j] в локальную префиксную сумму, так как правильные потомки - это PE с более высокими индексами
        send x[j] ⨁ m[j] to parent
    else:
        send x[j] to parent

// Нисходящая фаза
for j=0 to k-1:
    m[j] @ this = 0

    if hasParent:
        blocking receive m[j] @ parent
        // Для левого ребенка m[j] - родительская эксклюзивная префиксная сумма, для правого ребенка - инклюзивная префиксная сумма
        x[j] = m[j] ⨁ x[j] 

    send m[j] to left  // Общая префиксная сумма всех PE, меньших, чем этот или любой PE в левом поддереве
    send x[j] to right // Общая префиксная сумма всех PE меньше или равных этому PE

Конвейерная передача может быть применена в случае, когда сообщение ${\displaystyle m}$ длины ${\displaystyle n}$ может быть разделено на ${\displaystyle k}$ частей и оператор ⨁ может быть применен к каждой такой части отдельно.[16]

Если алгоритм используется без конвейерной передачи, то в дереве в любой момент времени работают только два уровня(отправляющий PE и принимающий PE), в то время как остальные РЕ ожидают. Если используется бинарное сбалансированное дерево ${\displaystyle p}$ обрабатывающих элементов, содержащее ${\displaystyle \log _{2}p}$ уровней, то длина пути от ${\displaystyle PE_{0}}$ до ${\displaystyle PE_{\mathbb {root} }}$ ${\displaystyle \log _{2}p-1}$, что соответствует максимальному числу непараллельных операций связи восходящей фазы. Аналогично связи по нисходящему пути также ограничены этим же значением. Считая время запуска связи ${\displaystyle T_{\mathbb {start} }}$ и байтовое время передачи ${\displaystyle T_{\mathbb {byte} }}$ получим, что фазы ограничены по времени ${\displaystyle (2\log _{2}p-2)(T_{\mathbb {start} }+n\cdot T_{\mathbb {byte} })}$ в неконвейерной передаче. При делении на ${\displaystyle k}$ частей, каждый из которых имеет ${\displaystyle {\frac {n}{k}}}$ элементов и отправляет их независимо друг от друга, первая часть будет требовать ${\displaystyle (\log _{2}p-1)(T_{\mathbb {start} }+{\frac {n}{k}}\cdot T_{\mathbb {byte} })}$ для прохождения к ${\displaystyle PE_{\mathbb {root} }}$ как часть локальной префиксной суммы и такое время будет действительно для последней части, если ${\displaystyle k>\log _{2}p}$.

В остальных частях все РЕ могут работать параллельно и каждая третья операция взаимодействия(получение слева, получения справа, отправка родителю) отправляет пакет на следующий уровень, поэтому одна фаза может быть проделана за ${\displaystyle 2\log _{2}p-1+3(k-1)}$ операций взаимодействия и обе фазы вместе требуют ${\displaystyle (4\cdot \log _{2}p-2+6(k-1))(T_{\mathbb {start} }+{\frac {n}{k}}\cdot T_{\mathbb {byte} })}$, что является очень хорошим показателем для сообщений длины ${\displaystyle n}$.

Алгоритм может быть дополнительно оптимизирован путем использования полнодуплексной или телекоммуникационной модели связи и перекрытия восходящей и нисходящей фазы.[16]