В работе рассматривается параллельная MPI+threads+SIMD реализация алгоритма вычисления разреженного QR разложения специальным образом упорядоченной прямоугольной матрицы на основе разреженных блочных преобразований лХ- аусхо дера. В алгоритме производится дпварреительное независимое параллельное высч-и ление QR разложений для наборов строк матрицы. Затем в соответствии с деревом зависимостей производится вычисление QR разложения матриц, составленных из окаймлений R- факторов строчных разложений. Приводятся рьетзаутлы эксперимнетов, подтверждающие эффективность предложенной параллельной реализации для тестовых задач. Алгоритм также может быть эффективно реализован нан-гетероге ных кластерных архитектурах с ускорителями типа GPGPU.
жении алгоритма на параллельную архитектуру компьютера. В Разделе 5 приводится описание
тестовой задачии представленырезультаты численных ксэпериментов.
В этом разделе приводится краткое описание параллельного алгоритма из работы [
Последовательный алгоритмпостроения QR разложениояснован на блочном преобрааз-ов нии Хаусхолдера вида , (2.1) где , , , . Блочное преобразование (2.1) строится через известный набор из обычных векторных преобразований Хаусхолдера следующим образом: матрица составляется из набора векторов направлений векторных преобразований Хл-аусхо дера, а верхняя треугольная матрицмаожет вычислена рекуррентным образом с испоал-ьзов нием коэффициенотв векторных преобразований Хаусхолдера если известна матрица .
При обсуждении разреженных вычислений будут рассматриваться вопросы вычисления QR разложения для прямоугольных так называемых мелко блочно разреженных матриц. Это означает, что разреежнность понимается в смысле блоков малого размера, каждый из которых является плотной в общем случае прямоугольной матрицей. Для простоты будем предполагать, что все мелкие блоки квадратные малого размер.а При этовмсе алгоритмы могут быть бо-бо щены на случай переменного столбцевого и строчного мелко блочного биения. В противовес мелким плотным блокам будем говорить также о кр-убплноочном или просто блочном биении, строчном и столбцевом. Это будет означать, что соответствующие подмосаттарви-цы с лены из некоторого количества мелко блочных строк и столбцов. При этом под е-блочным пр образованием Хаусхолдера имеется в виду преобразование вида (2.1) длямелокдонобглоочного столбца.
При последовательном построении профильного разреженногоQR разложения мелко блочной разреженной матрицы действуем по аналогии со случаем плотной матрицы.По превому мелко блочному столбцу матрицы строим разреженное блочное преобразао-в ние Хаусхолдера с разреженностью столбца такой, чтобы обнулитьл-ме кие блоки матрицы под первой блочной диагональю. Применяем транспонированное блочное прбеоразование Хаусхолдера ко второму мелко блочному столбцу, для уп-ол ченного результата строим слюед-у щее разреженное блочное преобазрова- Рис. 1. Профильное (слева) и расширенное профильное ние Хаусхолдера для обнуления элне-ме(справа) QR разложения тов под второй мелко блочноайгондаилью, и т.д.
Наравне с рпофильным разреженнымQR разложением рассмотрим также расширенное
профильное QR разложенимеатрицы. Схематически профильное и сшраиренное профильное
QR изображены на Рисунк1е. Фактически расширенное профильноQеR разложение - это
профильное QR разложение, примененное к матрице, расширенной сверху нулевым квадратным
блоком. При этом структура разреженноQстфиактора пополняется вмозожными
дополнитеьлными элементами на месте бывшего факRт,ориа отсоединенными диагональными элементами,
примыкающими к новомRу фактору.
В работе [
Для прямоугольнойразреженной матрицы введем ее строчное биение в виде: где , , и . Пусть для каждой из матриц имеет местоз-QR ра ложение с блочными преобразованиями Хаусхолд:ера
. Для разреженной матрицы ее строчные блоки могут содержатьмного нулевых
мелко блочных столбц,овв подобных случаях матрицы в 2(.3) не обязательно
верихентреугольные и имеют много нулевых стол,бцаов среди блочных преобразований Хаусхолдера
имеется многотождественных преобразований с единичной матрицей [
Рассмотрим задачу построения QR разложения всей матриц2.ы2) инза ( основе разложений (2.3). Для этого рассмотрим разреженную матрицу и ее QR разложение Обозначим . Имеют место соотношения (2.2) (2.3) (2.4) (2.5) (2.6) (2.7) ,
. Отсюда следует, что неявное представле2н.7и)е с(овместно с равенством из 2(.6) есть QR разложение матриц, ыпри этом матрицапредставляет собой -Q фактор QR разложения, а квадратная верхняя треугольная матерситцьа R- фактор QR раозлжения.
Описанная конструкция очевидно позволяет параллельным образом вычислять Qо-R разл жение матрицы за счет введения строчного б,иепноисякольук строчные QR разложени(я2.3) можно считать независимо. Синхронизация вычислений происходит только при вычислении объединяющего QR разложения2.5(). С другой стороны понятно, что подобный подсх-од к о новному распараллеливанию вычислений может быть эффенктитволько если число столбцов в матрице существенно меньше числа строк, иначе затраты на объединяющее QR разложение могут доминировать в вычислениях.
Описанный параллельный алгоритм вычисления QR разложения по блочным стржо-кам мо
но сделать эффективным сзчает использования дополнительной столбцевой разреженности
матрицы. Для этого рассмотрим д-вуурховневую организацию вычислений для прямоугольной
матрицы со структурой разреженности, изображенной на Рисунке 3. Пусть число мелко
блочных столбцов в мацтраих , и равны
соответтсвенно , и соответственно, . Как
показано в [
Определение 1. Прямоугольную мелко блочную миа-тр цу с введенным на ней блочным строчным и столбцевым биениями будем называтвьерхней квази- треугольной уровневой матрицей со структурой разреженности Рис. 3. Двух- уровневая органиазтипа вложенных сечений, если в терминах крупныох- блция параллельного вычислениQяR ков матрица является квадратной верхней треугольной и имеет структуру разреженности типа вложенных сечений, описывае-муорйовневым деревом зависимостей вычислений.
В частности, матрица на Рису3нкве терминах Определения 1 является двух уровневой верхней квази-треугольной с двух уровневым бинарным деревом зависимостей вычислений.
Как следует из предыдущего изложения, для ктэифвфноего вычисления QR разложения верхних квази-треугольных матриц с разреженностью типа вложенных сечений можньо- испол зовать следующий параллельный алгоритм: 1. Параллельно и независимо для каждой блочной строки матрицы строимн-расшире ное профильное разреженное QR разложение на основе блочных преобразований Хаусхолдера; 2. Параллельно в порядке, определяемом деревом зависимостей вычислений, и-достра ваем QR разложения для объединяющих подматриц для вычисления QRе- разлож ний соответствующих млкео блочных столбцевых окаймлений. 3. Архитектура гибридной вычислительной системы
Большинство современных суперкомпьютерных вычислительных систем как правие-ло им ют неоднородную архитектуру. С одной стороны, имеется набор вычислительных су-злов с ра пределенной памятью, обмен данными между которыми может быть осуществлен по быстрой обменной сети. С другой стороны, каждый узел представляет собой многро-процессо ный/многоядерный компьютер с общим доступом к оперативной памяти. Программин-ая реал зация вычислительных алгоритмов(включая алгоритм вычисления разреженного QR разел-ож ния) на компьютерах подобной архитектуры предполагает использование стандарта MPI при распараллеливании по распределенной памяти (по узлам вычислений), а по общей памяти узла распараллеливание по процессорам/ядрам естественно осуществлять на основе стандаа-ртов р боты с потоками с встроенными механизмами динамической балансировки нагрузки, таких как OpenMP или Intel® Threading Building Blocks (TBB). При этом предполагается4)(,Ричстуонок 4. Отображение алгоритма разреженного QR разложения на архитектуру вычислительной системы
Приведенный во втором разделе параллельнылгйориатм вычисления разреженного QR разложения для верхней кв-азтиреугольной матрицы типа вложенных сечений был реализован на кластерной MPI+threads архитектуре с использованием SIMD инструкций. Распаар-аллелив ние алгоритмана гетерогенной MPI+threads+SIMD архиттуерке осуществлено следующим боразом.
Распараллеливание верхнего уровня по MPI осуществлялось как распараллеливание по распределенной памяти. Для этих целей в дереве зависимостей вычислений каждомо-у MPI пр цессу было выделенцоеликом под- дерево зависимых вчыислений по возможности с близкой для всех поддеревьевычислительной работой. Дополнительная динамическая балансировка вычислений в како-лйибо форме не проводилась. Обрабоктакжадого извышестоящих узлов дерева зависимостей передавалась одномнуапр(имер первому) из тех процессоров, который обрабатывал один из узлов сыновдеайнного узла. На каждый MPI процесс перераспряе-дел лись те блочные строки матрицы, которые нужны для окончательной обработки своих узлов поддеревьев зависимостей вычислений.
Распараллеливание среднего уровня по нитям осуществлялось с помощью технологии Intel TBB либо как независимые, ликбаок зависимые вычисления. Зависимости описываются в виде под- графа зависимых вычислений для узлов поддеревьев своего MPI пр,онцеезсасвоирсаимые вычисления проводились при начальном вычислении QR разложенийблодчлняых строк. При проведении вычислений с узлдоемрева зависимых вычислений, не входящиMмPIв- поддеревья, вычисление объединяющих QR разложений проводилось только при поступлениио- необх димых днаных с других MPI процессов.
Распараллеливание нижнего уровня параллельных вычислен-иSйIMD векторизацияпроводилось за счет использования блочных преобразований Хаусхолдера. Рассмотрим этот вопрос подробнее. Основными операциями при работе с блочнпыремоибразованиями Хаусхолдера являются: 1. Вычисление векторного QR разложения с векторными преобразованиями Хаусхолдера для мелко блочного столбца;
2. Преобразование набора векторных преобразований Хаусхолдера в единое блео-чное пр образование Хаусхолдера дял мелко блочного столбца;
3. Применение с учетом разреженности блочных преобразований Хаусхолдерае- к посл дующим мелко блочным столбцам матрицы.
С учетом сказанного в Разделе 2спопсрооб трансформациивекторных преобразований Хаусхолдера в блочоне можно выделить следующие 4 основные вычислительные операции: векторный do:t - скалярное произведение векторов; векторный axpy: ; блочный do:t - блочное обобщение скалярногпороизведения для блоков ; блочный axp:y .
Векторные операции встречаются при векторном вычислении QR разло.жВенопиеярациях dot и axpдyлины векторов недостаточны для покрытия накладных расходов вызова ио-птимиз рованных BLAS функций из MKL. По этой причине в этих операциях осуществлялась ручная SIMD векторизация прямым вызовом соответствщуюих векторных инструкцисй помощью нитринсик функций.
Для максимальной локализации работы с памятью при обработке блочных преобразований Хаусхолдера естественно использовать формат хранения данных "по строкам" вмесцтио- тради онно используемого для блока векторов формата "по столбцам".
При проведении SIMD векторизацдилия блоков векторовв силу особенностей векторных инструкций реализовывалась поддержка только значений для двух типов данных float и double. Как уожтмеечалось, разреженные QR разложения будут использоватьсян- в ко тексте построения разреженных базисов в алгоритмах решения СЛАУ и задач наименьших квадратов, а в этом случае парам-етэрто выбранный размер блокиатерационной схемы.Для длинных блоков еквторов задачавычисления блочного dot и блочного axpy сводилась к циклу вызовов для подматриц размера . Подробности реализации операций блочного dot и axpy в терминах SIMD инструкций для подматриц стандартного размера можно найти в работе [6]. Как показали численные эксперименты, ручная векторизация для таких маленьких размеров блока оказалась значительно эффективней вызовов библиотечных реализаций из MKL и IPP. 5. Результаты численных экспериментов
Для тестирования предложенных в работе алгорвитмбыол выбран искусственный тест, в котором по возможности отражеонсыновные особенности будущегоих использования.
Для регулярной прямоугольной сетки была построернеагулярная разрежненая мелко блочная матрица с размерами бло,ков , сртуктура которой отвечает шаблону уравнения Лаплас.а Выбор размера блобкыал обусловлен тем, что для такого размера жв-озмо но для данных типов float и doпuрbоlдeемонстрировать возможности SIMD векторных т-инс рукций вплоть до набора AVX2.
Для полученной матрциы с помощью алгоритма вложенных сеч,екнаикй описано в работе
[
,
, arch\threads no-vec SSE AVX AVX2
arch\threads no-vec SSE AVX также как и, представляет собой мелко блочную с размером бл-оукраовневую верхнюю квази- треугольную матрицу типа вложенных сечений, и число столбцов евстьнепйолное чс-и ло столбцов бвазисе .
Эксперименты по вычислению разреженного QR разложения проводилиспьостдрлояенной таким образом матриц ы. В частности, если число векторов в бмазниосгео меньше мл-е ко- блочного размера матрицы, то в основных подматрицах число столбцов много меньше числа строк, что оправдывает использование расширенного разрежеQнRногроазложения вс- о новных вычислениях. Кроме того, каждый блобказиса обладает разреженностью, а потому результат произведения также есть разреженная матрица.
Для тестирования алгоритма была сгенерирована матрица с числом строк равным порядка 1,5 млн.и размером мелкого блока и количеством столбцевых блоков на два порядка меньше количества строк в матриТцеес.тирование алгоритма проводилось на-яд1е8рном процессоре Intel Xeon E5-2699v3 с архитектуройHaswell под управлениемCentOS 6.6. Тестовое приложение компилировалось с помощьоюптимизирующего компилятора ICC-15.0.3. В таблицах 1 и 2 представлены времена работы алгоритма для матрицы с одинарной и о-двойной точн стью соответственно для различных наборов векторных инструкций еиствкаолипчотоков. Таблица 1 Времена работы алгоритма с числами одинарной точности
Рис. 5. Профиль загрузки потоков в тестовом приложении. независимых блоков, а втор–аяокончанию обработки зависимых блоков по дереву. Видно, что при обработке зависимых блоков вычислений все еще остается значйитедлиьснбыаланс а-з грузки ядер, что не позволяет вплотную приблизкитьлсиянейной масштабируемости. В ь-дал Таблица 2 Времена работы алгоритма с числами двойной точности AVX2 1,043 0,528 0,282 0,166 0,12 0,0962 0,0904 Из результатов видно, что использование самых современных векторных инсптроу-кций зволяет получить ускорение до 3 раз по сравнению с оптимизирующим комIпCиCля.тоУср-ом корение по сравнению с бесплатным компилятGорCоCм получается еще более значительным. Использование всех 18 ядер процессора ускоряет работу алгоритма в обоих случаях в более чем 11 раз.На Рисунке 5 изображен профиль загрузки потоков в тестовом приложении, полученном с помощью программIыntel VTune Amplifier, оранжевым цветом отмечены региосинныхронизации потоков. Первая оранжевая область на временной линии отвечает окончанию обработки 4 0,609 0,312 0,245 0,196 4 0,677 0,555 0,399 нейшем планируется уделить больше внимания этим местам в алгТоаркиитмме.образом, с-и пользование обоих механизмов распараллеливания в современных проацхесспоорзволяет дс-о тичь ускорения в более чем 30 раз. Заключение
В работе представлаенреализация нагибридной параллельной MPI+threads+SIMD иар-х тектуре параллельного алгоритма вычисления QR разложенимяногоуровневой разреженной верхней квази- треугольной матрицысо структурой разреженносттиипа вложенных сечений. Результаты численных экспериментов с предложенным алгоритмом для тестовых гзиадба-ч на ридной параллельной MPI+threads+SIMD архитектурпеоказывают высокую эффективность предложенных алгоритмов: ускорение до 3 раз от использования векторных инсAтрVуXкц2и,й ускорение до 11 раз при использовании 18 ядер процессора. Литература 6. Применение векторных инструкций в алгоритмахчнбылхо операций линейной алгебры / Андреев А.Е., Егунов В.А., Насонов А.А., Новокщенов А.А. // Известия ВолгГТУ. Серия "Актуальные проблемы управления, вычислительной техники и информатики ве- технич ских системах".Вып. 21 : межвуз. сб. науч. ст. / В-оВлгоГлТгоУг.рад, 2014.- № 12 (139-). C. 5-11. Parallel implementation of the sparse QR decomposition for rectangular upper quasi triangular matrix with ND-type sparsity Sergey Kharchenko and Alexey Yushchenko Keywords: sparse rectangular matrix, upper quasi triangular matrix, nested dissection, QR decomposition, Householder transformations, MPI, multithreading, SIMD The paper considers parallel MPI+threads+SIMD implementation of the algorithm for computing sparse QR decomposition of a specially ordered rectangular matrix. Decomposition is based on block sparse Householder transformations. The algorithm starts with independent parallel QR decompositions for sets of matrix rows; and then, according to the computations tree, the QR decomposition is performed for matrices, combined with elements of R factors of rows decompositions. The results of numerical experiments for test problems show efficiency of the parallel implementation. The algorithm can also be efficiently implemented on heterogeneous cluster architectures with GPGPU accelerators.