Основная цель данной работы - разработка базовых вычислительных ядер, позволяющих перенести часть вычислительных функций программного комплекса Flow Vision [1] с CPU на GPU. Для этого необходимо было реализовать некоторый набор базовых операций линейной алгебры плотных и разреженных матриц, с учетом специфики данных, используемых в пакете FlowVision. Основной причиной невозможности использования уже готовых библиотечных функций от компании NVidia [2][3] было то, что они не предоставляют достаточно эффективных алгоритмов параллельного вычисления нескольких задач небольшой размерности на одной видеокарте. В результате данной работы было получено ускорение от 2 до 20 раз по сравнению с уже упомянутыми библиотечными аналогами.
Разработка базовых GPU ядер для новых блочных AMG алгоритмов для решения СЛАУ с явно вычисляемым базисом 3. Умножение блочной разреженной матрицы и транспонированной к ней на блок из плотных векторов Операция умножения блочной разреженной матрицы (или транспонированной к ней) на блок из плотных векторов представляет собой одну из двух следующих операций, то есть обновления блока векторов Y со знаком «+» или «-».
или (1.3) A – разреженная мелко-блочная матрица. Также умножение может производиться с учетом диагональных блоков матрицы или без. 4. (Не) полное треугольное разложение блочной разреженной матрицы Операция неполного треугольного разложения блочной разреженной матрицы представляет собой нахождение матриц, удовлетворяющих следующему отношению: или (1.4) где A, L, U – блочные разреженные матрицы с совпадающими профилями разреженности. 5. Решение СЛАУ Решение системы линейных уравнений с блочной разреженной верхней или нижней треугольной матрицей.
(1.5) Где A — блочная разреженная матрица, имеющая треугольную форму, B, X — блоки плотных векторов.
Все приведенные операции обладают рядом особенностей. Размеры блока матрицы и блока векторов принимают одно из дискретного набора значений: {4, 8, 16}. Они могут принимать как совпадающие, так и различные значения, однако первый случай наиболее интересен (например, операции DOT и DAXPY будут использоваться в операции QR разложения только с совпадающими размерами). Тип данных с плавающей точкой может быть как float, так и double, типы данных аргументов одной операции всегда совпадает.
Для эффективной реализации приведенных выше операций важно учитывать специфику хранения данных в пакете FlowVision. Во-первых, существуют приведенные выше ограничения на размер мелкого блока матрицы, блока из векторов, и.т.д. Благодаря этому можно проводить оптимизации, связанные с фиксированным размером блока (реализация отдельных ядер для каждого фиксированного размера). Во-вторых, все матрицы и вектора, используемые в вычислениях пакета FlowVision, обычно имеют небольшую размерность, вследствие чего невыгодно использовать всю видеокарту для вычисления одной операции. Решение этой проблемы как раз составляет главную идею представленной работы и описано в следующем разделе. 2. Общие результаты и исследования
Из-за особенностей организации вычислений пакета FlowVision (все вектора и матрицы имеют небольшие размеры) появляется проблема эффективного использования видеокарты, так как последовательно запускать множество задач небольшой размерности, занимающих всю видеокарту, как это делается при классическом использовании видеокарты, крайне невыгодно. В этом разделе описан основной подход, применявшийся для решения данной проблемы, а именно: выполнение каждой операции лишь на части ядер видеокарты, а не на всех. В результате этого подхода получена возможность эффективно запускать большое число параллельно работающих независимых операций на одной видеокарте, что и требуется от данных реализаций. Кроме того, при данном подходе появляется возможность эффективно использовать CPU и GPU одновременно.
Для эффективного запуска множества параллельных задач небольшой размерности на одной видеокарте необходима архитектура Kepler и выше, так как данные ускорители поддерживают технологию Hyper-Q. Данная технология необходима для того, чтобы можно было независимо запускать до 32 конкурирующих ядер в различных CUDA-потоках, что не было возможным в более ранних архитектурах. Каждое ядро представляет собой одну из 6 приведенных выше операций. Запущенное ядро использует GPU следующем образом: оно запускается в конфигурации 1 блок, 1024 нити. Благодаря этому, каждая операция вычисляется в одном из потоков непрерывно, после чего за ней может последовать другая операция (возможно, также другого типа). Далее последует объяснение, как данные ядра размещаются на мультипроцессорах видеокарты.
Благодаря тому, что каждая из операций занимает 1024 нити, на каждом мультипроцессоре могут одновременно выполняться одна или две операции (для архитектуры Kepler доступно максимально 2048 нити на мультипроцессор, причем максимальное число нитей в блоке 1024). На современных вычислительных GPU, используемых в суперкомпьютерных вычислениях, устанавливаются от 13 до 16 мультипроцессоров. Таким образом, при запуске двух задач на мультипроцессор на одной видеокарте можно запускать от 26 до 32 (в зависимости от модели видеокарты) параллельно работающих ядер. Благодаря технологии Hyper-Q, все запущенные таким образом ядра работают параллельно.
Исследование показало, что эффективнее размещать на каждом мультипроцессоре по две операции, так как загрузка вычислительных ядер мультипроцессора при таком подходе выше, что частично компенсирует потери из-за латентности доступа к памяти. Это подтверждается следующими диаграммами, демонстрирующими графики зависимости времени выполнения операций от числа параллельно запущенных на видеокарте задач (запускаются однотипные операции).
float double float double 1) 2) 20 я и16 н е н л12 о п в 8 ы я м е 4 р В
0 float double 1 2 4 6 8 101214151618202224262830
Число задач float
double 3) 4) 1 2 4 6 8 101214151618202224262830 Число задач Рис. 2.1. Зависимость времени выполнения(ms) от числа задач при фиксированном размере задачи для операций 1) DAXPY, 2) DOT, 3) и 4) матрично-векторного умножения.
Данные графики приведены для видеокарты K40, оборудованной 15 мультипроцессорами. Таким образом, при описанном подходе на ней можно запустить от 1 до 30 задач, причем при запуске 16-ой задачи время выполнения увеличивается меньше, чем в два раза. Очевидно, это происходит из-за того, что на одном мультипроцессоре запускается уже не одна, а две задачи. Кроме того, данный подход решает проблему одновременного использования CPU и GPU. Действительно, вычисления можно производить по следующей схеме: запускается некоторое количество CPU потоков (например, с помощью Intel TBB), некоторые из которых будут прикрепляться к ядрам CPU и проводить вычисления на них, в то время как другие будут заниматься управлением вычислениями на GPU (в основном запуском ядер и копированием данных). Эту схему демонстрирует следующая иллюстрация:
Рис. 2.2. Взаимодействие CPU и GPU задач через Intel TBB Таким образом, например, для 4-ядерного центрального процессора и видеокарты K40 с 15 мультипроцессорами можно запустить 3 + 30 = 33 CPU потока. Три из них будут вычисляться на ядрах центрального процессора (одно оставшееся свободное ядро будет выполнять управления видеокартой), остальные 30 на видеокарте. 3. Результаты по операциям 3.1. Операция DAXPY (3.1) Операция DAXPY, как уже говорилось, представляет собой прямое обобщение точечной операции AXPY. Где A – прямоугольный блок размера и M x , где M >> .
x , , а векторы x и y размера M x Операция реализована несколькими CUDA-ядрами для различных размеров мелкого блока. При совпадающих размерах получена наибольшая производительности из-за наиболее эффективного использования регистров и разделяемой памяти, а так же высокой занятости ядер графического процессора. Так же следует заметить, что наиболее интересны с практической точки зрения как раз случаи с квадратным мелким блоком размера 4x4, 8x8, 16x16.
Само ядро устроено следующим образом: для вычислений всегда используются 1024 нити, таким образом, при числе мультипроцессоров N можно параллельно запустить 2 * N таких ядер, как и описывалось в предыдущем разделе. Наиболее производительные ядра используют следующую стратегию вычислений: мелкий блок целиком помещается на регистры, а блок из векторов выгружается нитями из глобальной памяти с помощью двойного кэширования: сначала данные помещаются в разделяемую память, а затем на регистры. Умножение происходи только с использованием регистровых переменных. Так же, где возможно, используются операции __shfl() из архитектуры вычислений 3.0 для того, чтобы обойти разделяемую память и хранить данные только на регистрах.
Реализованные ядра можно сравнивать с библиотечной функцией DGEMM из cuBLAS, которая выполняет ту же операцию при правильном наборе параметров. Результаты сравнения их производительностей представлены на следующих диаграммах: cuBlas sgemm daxpy kernel cuBlas dgemm daxpy kernel 500 ь т со400 н ь ел300 т и од200 в з ои100 р П 0 16x16 4x4 8x8 Размер мелкого блока 16x16 Рис. 3.1. Сравнение производительности(GFlops) cuBlas и ядра DAXPY. 30 задач, симости от размера блока, точность float и double. элементов в завиИз приведенных выше диаграмм видно, что реализованные ядра DAXPY обгоняют библиотечные реализации от 4 до 20 раз. Достигается, это, в основном за счет параллельного использования видеокарты для вычислений, а так же оптимизаций, связанных с фиксированным размером блока. 3.2. Операция DOT
Реализует скалярное произведение блоков векторов. где А и B — блоки векторов, то есть матрицы размера МхN, где M > > N, а N — размер блока.
Результат С — блок размера NxN, такой, что C[i][j] — скалярное произведение i-го столбца матрицы А на j-й столбец матрицы B. В рамках исследования ядро запускается только на одном мультипроцессоре, занимая 1024 виртуальные нити, принадлежащие одному блоку. Так же, как и в операции DAXPY, используется двойное кэширование, сначала на быструю разделяемую память, потом на регистры. Все вычисления происходят на регистрах. Сравнение полученных результатов с операцией DGEMM из cuBLAS: cuBlas sgemm
dot kernel 250 ь т со200 н ьел150 т оид100 в зи 50 о рП 0 4x4 8x8 Размер мелкого блока 16x16 Рис. 4.1. Сравнение производительности(GFlops) cuBlas и ядра DOT. 30 задач, мости от размера блока, точность float.
элементов в зависи3.3. Блочное умножение на мелко блочную разреженную матрицу и транспонированную к ней
Умножение блочной разреженной матрицы на блок из плотных векторов (далее сокращенно MVM от англ. Matrix Vector Multiplication) представляет из себя одну из следующих двух операций, то есть обновление вектора Y со знаком «+» или «-» : (3.3) 1) Сначала к матрице добавляются две дополнительные структуры данных: массив структур, описывающий группы строк с одинаковыми длинами, и ссылки на новые индексы строк после перестановки. 2) Затем происходит перестановка индексов строк сортировкой длин с помощью стандартной функции quick sort на центральном процессоре. Затем, при необходимости вычислений на графическом ускорителе, все данные матрицы вместе с нововведенными структурами данных переносятся в память видеокарты. 3) После этого вычисления можно наиболее эффективно производить на графическом ускорителе.
Часть, выполняемая графическим процессором, состоит из следующих действий: 1) На графическом ускорителе запускается цикл по числу групп строк различной длины. 2) Для каждой группы строк вызывается __device__ функция, выполняющие умножение этих строк на вектор X. Для каждой строки выделяется группа нитей, с числом равным квадрату размера блока, однако, как уже говорилось – не более 1024 нитей на группу строк. 3) Блоки из каждой строки, умножаемые на соответствующие блоки из группы векторов, последовательно загружаются в разделяемую память и перемножаются. Благодаря удачно выбранному порядку циклов в умножении удалось использовать 100% пропускной способности разделяемой памяти, поэтому на регистрах хранится только промежуточная сумма. В конце обработки строки эта сумма сохраняется в глобальную память и группа нитей начинает обрабатывать следующую строку.
Для вычисления результата операции MVMT (умножения транспонированной разреженной матрицы на блок векторов) выполняются те же действия, но матрица предварительно переводится в формат CSC и вместо описания строк хранит описание столбцов матрицы. Это позволяет использовать то же ядро с минимальными модификациями. Благодаря такому подходу удалось эффективно реализовать операцию MVM на GPU даже для небольших размеров входных матриц (а так же для сильно разреженных матриц). Далее приведено сравнение реализованных ядер с библиотечным аналогом – функцией cusparseSbsrmv из библиотеки cuSparse. Результаты сравнения их производительностей представлены на следующих диаграммах: cuSparse kernel mvm kernel cuSparse kernel mvm kernel
Из приведенных выше диаграмм видно, что реализованные ядра MVM обгоняют библиотечные реализации от 2 до 6 раз. 3.4. Неполное треугольное разложение блочной разреженной матрицы В данной операции для блочной разреженной матрицы A необходимо найти две такие матрицы L и U, что в случае неполного и полного треугольного разложения соответственно. Матрицы A, L, U устроены так же, как в операции матрично-векторного умножения. Операция треугольного разложения, как на графическом, так и на центральном процессоре состоит из трех частей: 1) Подготовка данных о разреженности матрицы, их дополнение для дальнейшей обработки 2) Копирование данных (разреженности и блоков) на видеокарту (если хотим произвести вычисления на GPU, а не на CPU) 3) Запуск вычислительного ядра (или вычислительной функции на центральном процессоре) Из приведенного списка видно, что первая часть для CPU и GPU одинаковая и уже была реализована разработчиками из компании «Тесис», поэтому в данной работе рассматриваться не будет. Так же важно заметить, что с вычислительной точки зрения она гораздо менее затратная, чем сами вычисления, и, кроме того, в данной части много ветвлений и операторов перехода и небольшие ресурсы параллелизма, поэтому выгоднее выполнять её на центральном процессоре.
Далее приведено описание схемы работы вычислительного ядра на графическом процессоре.
Вычислительное ядро запускается на 512 GPU-нитях (а не на 1024, как для предыдущих операций), так как оно использует значительно больше регистров, что ограничивает occupancy.
Для каждой строки матрицы выполняется: 1) Обновление текущей строки через предыдущие строки верхнего треугольника, а так же столбцы нижнего треугольника Данная функция вычисляется следующим образом: на каждую строку выделяется группа из всего множества нитей (так как нитей может быть сильно больше, чем необходимо для обработки одной строки при небольшом размере блока матрицы) – например 64 нити из 512 при обработке матрицы с блоком 4x4. Важно заметить, что эта часть вычислений в большинстве случаев наиболее затратная по времени выполнения. 2) Модификация диагональных блоков при необходимости (только в случае не полного разложения) Производится быстро через __shfl инструкции. 3) Факторизация и обращение диагонального блока Здесь диагональный блок, представляющий собой плотную матрицу размера блока (NxN), сначала приводится к треугольному виду, а затем полученный треугольник обращается (в симметричном случае). В обоих случаях параллельно работать может только N нитей, в то время как остальные вычисления производить нельзя, что сильно замедляет вычисления при небольшом числе блоков в строке. Полученный блок передается в следующую функцию через разделяемую память, что позволяет уменьшить число обращений к более медленной глобальной памяти. 4) Умножение блоков данной строки на обращенный блок Умножение всех блоков строки на обращенный блок из предыдущей операции. В обоих случаях сделано так, что каждая пара умножаемых блоков хранится на регистрах. 5) Сохранение результата в нижнюю и верхнюю треугольные матрицы Сохранение результатов в глобальную память (блоков верхнего и нижнего треугольных факторов). 3.5. Решение СЛАУ Данная операция представляет собой решение следующей системы: (3.5) где А – мелко-блочная разреженная треугольная матрица, В – плотный блок векторов, матрицы А хранятся по столбцам, а матрицы В – по строкам.
В этой задаче требуется вычислить не решение системы линейных уравнений с блочной матрицей А, а найти решения для блока систем линейных уравнений. Матрицы этих систем сгруппированы поэлементно в блочную матрицу А, а правые части сгруппированы в блок векторов В. Так как матрица уже дана в треугольной форме, то реализуется обратный ход метода Гаусса решения системы линейных уравнений, с той лишь разницей, что вместо элемента матрицы и элемента правой части используются блоки.
1. Сначала матрица переводится из формата CSR в формат CSC для удобства и локальности обращений в память.
2. Блок нитей работает с последним, не обновлённым столбцом, при этом параллельно вычисляются новые значения правых частей. Доступ к памяти осуществляется блоками, реализуется кэширование сначала в разделяемую память, оттуда на регистры. Все вычисления проводятся на регистрах. 4. Заключение
В ходе данной работы были эффективно реализованы шесть базовых операций линейной алгебры из пакета FlowVision. Для каждого из реализованных ядер было проведено большое количество оптимизаций, исследованы различные их характеристики и свойства. Также для каждой операции была теоретически обоснована эффективность реализации. Из анализа эффективности c использованием профилировщика Visual Profiler было определено, что основным ограничением на производительность написанных ядер оказалась пропускная способность разделяемой памяти, которая используется в качестве КЭШа с быстрым доступом для хранения обрабатываемых элементов матрицы (обычно мелкого блока). 1 200
GPU speedup
Dot MVM
Рис. 4.1. Ускорение относительно одноядерных тестовых CPU аналогов данных операций. Таким образом, в будущем эти ядра могут быть успешно использованы в программном пакете FlowVision компании «Тесис». Так же очень важно, что было получено ускорение в 2-20 раз относительно функций из программных пакетов cuBlas и cuSparse от компании NVidia, реализующих требуемые операции.
Работа выполнена/выполняется при финансовой поддержке Министерства образования и науки Российской Федерации, Соглашение №14.607.21.0006, уникальный идентификатор RFMEFI60714X0006. Литература The development of basic GPU kernels for the new block AMG algorithms for solving SLE with explicitly calculated sparse basis Ilya Afanasyev and Yury Potapov The main purpose of this work is the development of basic kernels, which make it possible to transfer some computing functions of software package Flow Vision from CPU to GPU. It was necessary to implement a set of basic dense and sparse linear algebra operations, with the data specifics used in the package FlowVision. We could not use already-created by company NVidia library functions because they do not provide sufficiently effective algorithms of computing several small dimension tasks on a single GPU. As a result of this work acceleration from 2 to 20 times was obtained in comparison with the already mentioned library counterparts.