Сортировка на GPU. Битоническая сортировка

Во многих случаях при рендеринге объектов (не говоря уже о GPGPU) возникает необходимость сортировки объектов (это может быть сортировка частиц в системе частиц, сортировка полупрозрачных граней и т.п.).

Однако традиционные методы сортировки к GPU неприменимы, а передавать данные на CPU для их сортировки крайне нежелательно.

Поэтому возникает необходимость в методах сортировки данных, хорошо ложащихся на архитектуру современных GPU.

Существенной особенностью таких алгоритмов является их многопроходность и параллельный характер. Т.е. при помощи фрагментного шейдера и рендеринга в текстуру строится последовательность текстур с данными. На последнем проходе мы получаем на выходе уже полностью отсортированный массив (текстуру).

Поэтому важным показателем для таких методов является число проходов (render passes), которые нужно выполнить для сортировки n элементов (далее мы будем считать что число элементов является степенью двух).

Рассмотрим сначала случай, когда у нас есть одномерный массив с данными (одномерная текстура), которую необходимо отсортировать.

Простейший вариант заключается в том, что на каждом проходе мы сравниваем между собой элементы x_2i и x_2i+1, переставляя их, если они не удовлетворяют принятому критерию упорядочения (далее мы будем считать, что мы хотим отсортировать массив по возрастанию).

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

Однако такому методу (называемому odd-even transition sort) для сортировки массива из n элементов потребуется ровно n проходов (чтобы убедится в этом достаточно рассмотреть случай, когда наименьший элемент находится в самом конце массива - за один проход он будет опускаться ровно на одну позицию).

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

И тут нам на помощь приходят понятия сетей компараторов (comparator networks) и сортирующих сетей, которые очень хорошо ложатся на архитектуру GPU.

Компаратор можно представить себе как устройство с n входами и n выходами, которое выполняет некоторую перестановку элементов в зависимости от результатов их сравнения.

Рис 1. Компаратор.

Из компараторов можно построить сеть, на вход которой поступает последовательность a₀,...,a_n-1, а на выходе получается последовательность b₀,...,b_n-1. Размером такой сети называется число компараторов в ней, а ее глубиной - максимальное число компараторов, проходимое одним входным элементом.

Рис 2. Пример компаратора.

Одной из простейших сетей компараторов является так называемый полуочиститель (half-cleaner) B_n, определяемый следующим образом: он сравнивает между собой элементы x_i и x_i+n/2 для всех i от 0 до n/2-1.

Рис 3. Компаратор B_n.

Т.е. сеть B_n сравнивает между собой элементы x_i и x_i+n/2 и упорядочивает их заданным образом (по возрастанию). При этом каждый элемент обрабатывается всего один раз (т.е. глубина сети равна единице).

На следующем рисунке приведены результаты применения к массиву (6,5,3,0,2,4,7,1) сетей B₈, B₄ и B₂.

Рис 4. Применение компараторов B₈, B₄ и B₂.

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

Обратите внимание, что сеть B_n очень легко ложится на архитектуру GPU и требует ровно одного прохода.

Можно путем применения определенной последовательности полуочистителей отсортировать произвольный массив заданного размера, т.е. из отдельных полуочистителей можно построить сортирующую сеть.

Справедливо следующее утверждение (его полное доказательство можно найти в Алгоритмы: построение и анализ):

Утверждение. Если сеть компараторов упорядочивает (сортирует) произвольную последовательность из нулей и единиц, то она является сортирующей, т.е. она упорядочивает произвольную последовательность чисел.

Это утверждение позволяет заметно упростить анализ сортирующих сетей, сведя все к наборам из нулей и единиц. Введем следующее определение.

Определение. Последовательность a₀,...,a_n-1 называется битонической, если она сперва убывает, а потом возрастает, либо получается из такой последовательности циклическим сдвигом.

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

Битоническая последовательность из нулей и единиц имеет вид 1,..,1,0,...,0.,1...,1 либо 0,...,0,1,..,1,0,...,0.

Справедливо следующее свойство битонической последовательности:

Утверждение. Если применить сеть B_n к битонической последовательности a₀,...,a_n-1 (длины n) из нулей и единиц, то получившееся в результате последовательность b₀,...,b_n-1 будет обладать следующими свойствами:

1. Ее верхняя и нижняя половины также будут битоническими.

2. Любой элемент первой (верхней) половины всегда будет меньше или равен любого элемента второй (нижней) половины.

3. Хотя бы одна из половин - чистая (т.е. монотонная).

Пусть у нас есть битоническая последовательность a₀,...,a_n-1. Если мы применим к ней сеть B_n, то мы в результате получим две половины, каждая из которых является битонической и все элементы первой половины меньше или равны элементов второй, т.е. если мы отсортируем каждую из половин по-отдельности, то мы получим отсортированную последовательность.

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

К каждой из этих четырех частей применим сеть B_n/4, в результате чего мы получим восемь правильно упорядоченных по отношению друг к другу частей, каждая из которых является монотонной.

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

Таким образом мы фактически построили сортирующую сеть для битонических последовательностей. Подобная сеть также называется битоническим слиянием (bitonic merge). Ее структура представлена на следующем рисунке:

Рис 5. Битоническое слияние.

Видно, что операция битонического слияния M_n может быть представлена рекурсивно в следующем виде (см. рис 5.).

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

Как легко видно из рисунков операция битонического слияния M_n требует ровно log₂n проходов.

Рассмотрим теперь как можно при помощи сетей B_n и M_n построить сеть, сортирующую произвольную последовательность длины n.

Пусть у нас есть последовательность a₀,a₁,a₂,a₃,a₄,a₅,a₆,a₇. Разобьем ее на пары и применим к каждой паре сеть B₂, но с чередующимся порядком сортировки.

В результате мы получим последовательность со следующим порядком (стрелки обозначают возрастание/убывание).

Рис 7. Результат применения B₂ с чередующимся порядком упорядочивания.

Как легко видно и первые четыре элемента (a₀,a₁,a₂,a₃) и вторые четыре элемента (a₄,a₅,a₆,a₇) получившейся последовательности являются битоническими. Поэтому если к первой применить операцию битонического слияния M₄, сортирующей по возрастанию, а ко второй - операцию слияния M₄, сортирующей по убыванию, то мы в результате получим последовательность, которая является битонической, т.е. состоящей из двух монотонных участков.

Рис 8.

Если теперь применить ко всей последовательности операцию M₈, то в результате мы получим полностью отсортированную последовательность.

Это последовательность операций и называется битонической сортировкой и ее рекурсивное представление изображено на следующем рисунке.

Рис 9. Рекурсивное представление битонической сортировки.

Как несложно убедиться, битоническая сортировка требует log₂n*(log₂n+1)/2 проходов.

Реализация битонической сортировки на GPU

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

Для простоты будем считать, что это массив 4-мерных векторов, которые мы будем сортировать по z-компоненте.

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

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

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

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

Если же элемент является вторым в паре, то в gl_FragColor записывается элемент с максимальным значением z-координаты.

Проще всего эту дополнительную информацию для шейдера передать в z и w компонентах текстурных координат.

Условимся, что отритцательные значения gl_TexCoord [0].z соответствуют первому элементу пары, а положительные - второму элементу.

Также условимся, что отритцательные значения gl_TexCoord [0].w соответствуют сортировке по убыванию, а положительные - сортировке по возрастанию.

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

В реализации управляющего кода (рендеринга) удобно произвести группировку отдельных проходов в стадии (stages). Каждая стадия фактически отвечает за слияние группы битонических последовательностей (с предыдущей стадии).

Так первая стадия будет состоять только из операции B₂ с чередующимся упорядочиванием. На выходе этой стадии мы получим набор пар с чередующимся способом упорядочивания (возрастание-убывание-возрастание- и т.д.).

Следующая стадия состоит из двух проходов - B₄ и B₂ - дает на выходе сортированные четверки элементов с чередующимся порядком упорядочивания (каждые две пары объединяются в одну четверку).

Из всех стадий выделяется только последняя (соответствующая операции M_n ) поскольку в ней всюду используется одно и тоже направление упорядочения. Поэтому ее удобно отделить от всех остальных стадий.

Как несложно убедиться всего у нас будет log₂n стадий.

Таким образом вся процедура битонической сортировки одномерного массива может быть представлена следующим образом:

Поскольку мы будем довольно много раз осуществлять рендеринг в текстуру (FBO), то удобно сразу же завести массив из двух FrameBuffer'ов и целочисленную переменную index, которая будет соответствовать тому буферу в массиве, который содержит данные для очередного прохода, а index^1 будем соответствовать буферу, куда очередной проход будет осуществлять рендеринг.

Рассмотрим сначала функцию doLastStage как более простую (по сравнению с функцией doStage, реализующей остальные проходы).

Как легко убедится в стадии в номером stage (stage принимает значения 1,...,log₂) будет ровно stage проходов.

Последний проход выполнят следующие операции - B_n/2^stage,...,B₄,B₂.

Один проход для последней стадии применяет операцию B_n к частям массива.

Рис 10. Пример последней стадии для битонической сортировки восьми элементов.

Как видно из приведенного рисунка на каждом проходе полуочиститель B_n применяется не ко всему массиву, а лишь к некоторому набору его элементов.

Для B₈ в качестве такого блока выступает весь исходный массив. Для B₄ отдельно обрабатывается первая четверка элементов и вторая четверка. Для полуочистителя B₂ отдельно обрабатываются четыре пары элементов.

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

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

Рассмотрим теперь как происходит обработка остальных стадий. На следующем рисунке приводится вторая стадия для сортировки массива из восьми элементов.

Рис 11. Вторая стадия для сортировки массива из восьми элементов.

Обратите внимание каким образом происходит смена направления сортировки - она фактически задается первым проходом стадии.

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

Так для второго прохода этой стадии будет два блока из четырех элементов (в первом сортируется по возрастанию, а во втором - по убыванию). И каждый такой блок разбивается на два подблока по два элемента.

Ниже приводится реализация данной функции.

Полностью текст программы (как и уже откомпилированную версию) можно скачать по ссылке в конце статьи.

Реализация битонической сортировки для двухмерных текстур

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

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

Гораздо перспективнее другой подход - мы сразу сортируем двухмерную текстуру. К сожалению первая приходящая в голову мысль - отсортировать матрицу по строкам, а потом по столбцам (благо код для этого очень легко получается из кода для одномерной сортировки) не работает (достаточно рассмотреть матрицу a_ij=i+j).

Понятно, что можно легко адаптировать код одномерной сортировки для сортировки текстуры по строкам (фактически нужно только поправить y-диапазон в коде рендеринга). Однако набор одинакового отсортированных строк нам мало что дает.

Гораздо интереснее будет, если порядок сортировки чередовать - первую строку сортировать по возрастанию, вторую - по убыванию, третью - опять по возрастанию и т.д.

Тогда каждая пара строк будет образовывать битоническую последовательность, которую можно объединить при помощи операции bitonic merge. На самом деле мы получим полный аналог сортировки одномерного массива после применение B₂ с чередованием порядка сортировки.

Рис 12. Строки с чередующимся порядком упорядочения элементов.

Только в этом случае мы будем уже объединять не группы элементов, а группы строк. Фактически это будет той же битонической сортировкой для очень большого одномерного массива, просто размеры объединяемые блоки будут представлять из себя наборы строк. В результате мы действительно получим отсортированный массив с той же оценкой числа проходов (только в качестве n будет выступать общее число элементов в матрице).

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

На следующих двух листингах приводится исходных код для этих шейдеров.

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

Общая процедура битонической сортировки выглядит следующим образом:

При этом функция rowSort осуществляет сортировку внутри строк с чередующимся порядком и выглядит следующим образом:

Как видно из приведенных листингов появились версии функций выполнения стадий для строк и столбцов, ниже приводятся реализации лишь двух из них.

По этой ссылке можно скачать весь исходный код к этой статье. Также доступны для скачивания откомпилированные версии для M$ Windows, Linux и Mac OS X.

Ниже приводится небольшой список литературы по этой теме.

1. A Memory Model for Scientific Algorithms on Grpahics Processors. Naga K. Govindaraju, Scott Larsen, Jim Gray, Dinesh Manocha.

2. A Cache-Efficient Sorting Algorithm for Database and Data Mining Computations using Graphics Processors. Naga K. Govindaraju, Nikunj Raghuvanshi, Michael Henson, David Tuft, Dinesh Manocha.

3. GPU-ABiSort: Optimal Parallel Sorting on Stream Architectures. Alexander Gres, Gabriel Zachmann.

4. GPU Gems II. Chapter 46. Improved GPU Sorting. Peter Kipfer, Rudiger Westermann.