CC BY
173
УДК 519.613.2
ВЫСОКОПРОИЗВОДИТЕЛЬНЫЙ АЛГОРИТМ ШЕРМАНА - МОРРИСОНА ОБРАЩЕНИЯ МАТРИЦ НА GPU 1
H.С. Недожогин, А.С. Сармакеева, С.П. Копысов
Обращение матрицы является важным этапом при численном решении таких, задач как решение систем линейных уравнений и построение предобуславливателей, вычисление дополнения Шура в методах декомпозиции области, цифровая обработка изображений и т. д. Разработка высопроизводительных параллельных алгоритмов обращения матриц связана с эффективным хранением и отображением алгоритмов на современные многоядерные архитектуры. Наряду с традиционными методами обращения — LU-факторизацией и методом Гаусса - Жордана, рассмотрены параллельные алгоритмы метода сопряженных градиентов и Шермана - Моррисона, в которых используются матрично-векторные и скалярные произведения эффективно выполняемые на многоядерных процессорах. В работе проведено сравнение на тестовых матрицах рассматриваемых методов на CPU и GPU.
Ключевые слова: высокопроизводительные алгоритмы, обращение матриц, разреженные матрицы, алгоритм Шермана - Моррисона.
Введение
При решении практических задач возникает необходимость выбора того или иного метода обращения матриц или его параллельной реализации. Появление программного обеспечения для вычислений общего назначения на графических устройствах (GPGPU) позволило на ряде задач, в том числе вычислительной линейной алгебры, получать ускорение вычислений в десятки и сотни раз, по сравнению с центральным процессором (CPU). Тысячи потоков GPU могут эффективно выполнять одновременно большое число простых арифметических операций, что характерно для мультипликативных и аддитивных операций с векторами и матрицами. Вместе с тем, последовательные операции и ветвления, характерные для разложения матриц на треугольные множители, выполняются GPU медленнее, чем ядрами CPU.
Данная работа ориентирована на эффективные реализации методов обращения симметричных разрежённых и положительно определенных матриц с использованием гибридных вычислительных узлов с GPGPU. В работе проведено сравнение алгоритмов обращения матриц, имеющих разреженную структуру. При реализации всех методов учитывалась структура хранения в сжатом формате.
I. Алгоритмы обращения матриц
Рассмотрим высокопроизводительные алгоритмы явного и неявного обращения матриц на гибридной архитектуре с выделением операций эффективно выполняемых на универсальных и графических процессорах. Хотя явное обращение матрицы во многих случаях может быть нецелесообразным, оно по-прежнему представляет интерес в различных прикладных областях, где встречаются большие разреженные и, порой, плохообусловленные матрицы.
1 Статья рекомендована к публикации программным комитетом научной конференции «Параллельные вычислительные технологии — 2014».
Ьи-факторизация (Ьи). Для заданной невырожденной матрицы А существует разложение на верхнюю треугольную и и на нижнюю треугольную Ь матрицы. Если матрицы А, и, Ь обратимы, то справедливо следующее утверждение: А-1 = и-1Ь-1. Так как факторизация проводилась для матриц хранящихся в формате СБИ, при построении Ь и и суммировались лишь те произведения, в которых участвовали ненулевые элементы исходной матрицы [1]. Рассматривалась только последовательная реализация Ьи-факторизации на СРи, показывающая наилучшие результаты для малых и средних матриц. Особенности данного алгоритма не позволяют эффективно использовать все возможности массивнопараллельной архитектуры ОРИ.
Метод Гаусса — Ж^ордана (О^. Потенциально, метод Гаусса - Жордана больше подходит для вычислений на графических ускорителях, чем алгоритм основанный на Ьи-разложении. В рассматриваемом варианте все операции выполняются на центральном процессоре СРи. Так как преобладающие вычисления в данном методе содержат матричные операции, то в дальнейшем можно ожидать высокую производительность варианта для гибридной архитектуры при совместном использовании СРИ+ОРИ. Однако, в работе [2] показано, что максимально получаемое трехкратное ускорение метода обращения Гаусса -Жордана достигается лишь при вычислениях с одинарной точностью и разделением операций выполняемых на СРИ и ОРИ.
Метод сопряженных градиентов (СС). Рассматриваемый в работе алгоритм вычисления обратной матрицы состоит из решений матричной системы вида АХ = Е, где Е
— единичная матрица. Такая система эффективно решается на ОРИ предобусловленным алгоритмом сопряженных градиентов. Использование в тестах простейшего диагонального предобуславливателя связано с обеспечением минимальных вычислительных затрат (количества требуемой памяти, времени построения предобуславливателя и времени его решения) на ОРИ для матриц общего вида.
Реализация алгоритма более подробно рассмотрена в работе [3]. Отметим только, что оценка числа ненулевых элементов обратной матрицы затруднена и обращенную матрицу приходится хранить полностью, что накладывает существенные ограничения на используемую память при работе с ОРИ. В случае обращения матрицы методом сопряженных градиентов результат вычислений получается по столбцам, что позволяет сразу преобразовывать матрицу в сжатый формат и, тем самым, снижать ограничение на размер задачи, решаемой на ОРИ. В методе сопряженных градиентов операции скалярного произведения, нормы векторов, суммы векторов и умножения вектора на скаляр распараллелены с помощью библиотеки сиБЬЛВ. Матрично-векторное произведение реализовано с помощью технологии СИБЛ в векторном варианте с возможностью использования нескольких ОРИ.
Метод Шермана — Моррисона (БМ). Идея метода БМ восходит к работам [4, 5], который известен в русскоязычной литературе метод как метод пополнения [6]. Сначала рассмотрим тождество Вудбери [4] из которого следует известная формула Шермана -Моррисона. За основу построения матрицы А-1 возьмем матрицу В той же размерности, что и А, но с известной обратной матрицей.
Теорема 1. (Вудбери [4]) Для обратимости матрицы А вида А = В — иУт, где А — невырожденная матрица порядка п х п, а и, У пусть (п х т) — матрицы, необходимо и достаточно, чтобы была обратимой матрица т-го порядка Р = 1т — УтВ-1 и. При этом
А-1 = В-1 + В-1иР-1УТ В-1. (1)
В частном случае, при т = 1, тождество Вудбери приводит к следующей теореме:
Теорема 2. (Шерман - Моррисон [5]) Пусть В — невырожденная матрица и вектора и и V такие, что
г = 1 + vT В-1и = 0.
Тогда матрица А = В-1 + vTи является обратимой и её обращение находится как
А-1 = В-1 - г-1 В-1 uvTВ-1. (2)
Далее будем рассматривать и считать, что матрица А обратимая с ац = 0 и нет необходимости в перестановках при её обращении. В этом случае, при условии и = вд =
[0... 1... 0]т, выражение (2) можно упростить, полагая что В-1вд = ЬД, ЬД — к-ый столбец
матрицы В-1. Тогда столбец а* матрицы А-1 вычисляется
vT Ь*
а* = Ь* - !+^ЬД’ * = 1 2’...’^. (3)
По соотношениям (3) вычислим обращение матрицы А, принимая что В = Е и записывая на к шаге
и = вд, VT = vk = ак - вт, (4)
исходя из
N
гк
A = в + £-
/ V
к=1
Столбец матрицы на к шаге находится из соотношения
^ а(к-1)
а(к) = а(к-1) і а(к-1) А = 12 N (5)
аі = аі “ , т (к-1) ак , А = Х ,2 , •••^, (5)
1 + Vі ак
где А(0) = Е и для к = N находится а(М) = а*, А = 1, 2,..., N.
Реализация алгоритма (2) в виде (5) приведена в Алгоритме 1. Векторные операции реализовывались на ОРОРИ с помощью функций библиотеки сиБЬЛВ.
Алгоритм 1 Параллельный Шермана - Моррисона
2
3
4
5
6
7
8 9
10
11
12
Require: a, v £ RN {вектора хранятся в памяти GPU}
1: vT = ak — eT { вектор efc £ Rn , k элемент которого равен 1, а все остальные 0} dim {размерность матрицы } for k = 0 ^ N do for i = 0 ^ N do for j = 0 ^ N do
в1 — (vT, ajk ^) {вычисляется с помощью функции cublasDdot} в2 — 1 + (vT, akk-1)) {cublasDdot}
в3 — в1/в2
a(k) — ajk-1) — в3 {cublasDaxpy} end for end for end for
В работе [7] был предложен блочный алгоритм обращения матриц на основе (1) в случае плотных матриц, в котором требуется достаточно затратное разделение матриц на блоки. Ускорение при использования GPU в этом случае составило лишь 20%.
В настоящей работе рассматривается алгоритм обращения матриц с использованием формул (2) для разряженных симметричных матриц.
2. Численные исследования
Численные эксперименты, представленные в таблице, были проведены на тестовых матрицах известных библиотек [8, 9], а также на матрицах полученных при решении задачи напряженно-деформированного состояния винтовой пружины методом дополнения Шура [3]. Незаполненные поля таблицы соответствуют вычислительным затратам существенно превышающим результаты по другим методам. Размеры матриц варьируются от 276 до 18000, число ненулевых элементов — от 1666 до 6897316, числа обусловленности %
— от 10 до 1010. Все результаты получены с двойной точностью.
Для вычислений использовался узел с четырехядерным процессором Intel Xeon CPU E5430 частотой 2667 МГц и графическим ускорителем NVIDIA GeForce GTX 580. Количество ядер CUDA 512, объем памяти 3 ГБ, частота ядра/памяти 772 МГц/4008 МГц.
Выделим три основных параметра для сравнения и оценки алгоритмов обращения матриц на параллельной вычислительной системе: число арифметических операций, затраты памяти и требования к обработке данных.
Затраты памяти метода сопряженных градиентов и метода пополнения на формирование исходной разреженной и обратной матриц составляют 16NNZ + 4N2 и 8N2 байт, соответственно. В методе Шермана - Моррисона на вспомогательные переменные, требуется 16N, в методе сопряженных градиентов затраты несколько больше и составляют 48N. Доступ к элементам матрицы в формате CSR осуществляется по столбцам в методе SM и по строкам в CG, тогда как в методах Гаусса - Жордана и LU-факторизации доступ к матрице происходит поэлементно и на каждом шаге необходимо иметь всю матрицу. Вычислительные затраты алгоритмов по числу операций выглядят следующим образом: в методе сопряженных градиентов без предобуславливания на каждой k-итерации один раз выполняется матрично-векторное произведение и три скалярных произведения N(N2 + 3N + 3); в методе пополнения на каждой итерации N раз выполняется два скалярных произведения N2(2N + 2).
Для выполнения скалярных произведений в методе SM создается вектор в памяти GPU, в который записывается столбец матрицы, необходимый на данном шаге алгоритма. В матрично-векторном произведении для метода CG каждую строку матрицы обрабатывает несколько нитей (до 32) одновременно, что обеспечивает его выполнение в два и, даже, три раза быстрее скалярного произведения при произвольных размерах матриц и их разрежённости.
Определяющей, при использовании метода сопряженных градиентов, становится обусловленность (х) обращаемой матрицы. В случае с хорошо обусловленными матрицами (когда число итераций метода сопряженных градиентов меньше размера матрицы), число операций в методе сопряженных градиентов будет меньше, чем в методе SM (см. табл. 1).
При плохой обусловленности матриц, но большой разреженности (Dubcova1 16129/253009, bodyy4 17546/121550) CG работает на порядок быстрее, чем метод SM. Од-
Таблица
Время обращения матриц, с
GJ LU CG SM
Название N/Nnz X CPU CPU GPU CPU GPU
Schur276 276/7488 1.2e+02 0.1 0.1 1.1 0.1 4.1
Lect 494_BUS 494/1666 3.9e+06 0.3 0.7 23.9 1.31 12.1
1138_BUS 1138/4054 8.0e+09 9.5 11.9 133.1 23.0 62.5
Schur1236 1236/40806 1.5e+01 12.6 18.6 8.8 34.5 74.5
BCSSTK11 1473/34241 5.3e+08 21.4 33.9 246.6 57.9 104.8
BCSSTK15 3948/117816 8.0e+09 426.0 649.1 329.1 1182.9 751.1
Shur5688 5688/203238 2.9e+01 1262.8 1822.5 84.8 3171.6 1572.2
ND3K 9000/3279690 1.6e+07 5009.6 7175.4 34801.3 12960.3 3918.7
msc10848 10848/1229776 9.9e+09 8718.8 12446.8 15259.6 21663.5 5741.0
Dubcova1 16129/253009 9.9e+02 29868.9 40859.7 617.3 75472.2 13006.5
bodyy4 17546/121550 8.1e+02 40400.6 — 376.4 — 15566.6
ND6K 18000/6897316 1.5e+07 41782.0 35674.5 — — 16386.0
нако, с увеличением числа ненулевых элементов затраты возрастают (ND3K 9000/3279690, ND6K 18000/6897316), а девятикратное ускорение достигается в параллельном методе SM.
В численных экспериментах, для выполнения базовых операций рассматриваемых методов, вызывались функции библиотеки cuBLAS (см. Алгоритм 1). Особенностью применения функции cublasDdot является то, что при использовании текстурной памяти в основной программе последующий вызов функции cublasDdot выполняется существеннее медленнее, чем другие вызовы при вычислении скалярных произведений. С увеличением размерности векторов это влияние становится значительным и операция скалярного произведения выполняется, практически, в семь раз медленнее, что характерно для метода CG, при реализации которого текстурная память задействована в матрично-векторном произведении.
Используемая реализация алгоритма CG хорошо масштабируется на GPU и позволяет достигать ускорения в десятки раз, по сравнению с вариантом, реализованном на CPU [1], а параллельная реализация алгоритма SM на GPU обеспечивает гораздо меньшие ускорения. Так, для матриц малых размерностей (Schur276, 494_BUS, 1138_BUS, Schur1236, BCSSTK11) время выполнения варианта, реализованного на GPU, больше, чем на CPU, а на матрицах больших размерностей получено максимальное шестикратное ускорение. Прежде всего это связано с лучшей масштабируемостью скалярного произведения небольших векторов на CPU, а не на GPU. В случае метода сопряженных градиентов значительное ускорение вычислений на GPU обеспечивает матрично-векторное произведение, которое обладает большой степенью параллелизма.
Заключение
При использовании GPU на малых задачах (до N = 3000) сокращение времени выполнения операций не покрывают затраты на выделение памяти, поэтому эффективность параллельных реализаций на CPU в два раза выше. Однако, с увеличением размеров матриц, алгоритмы реализованные на CPU не показывают значимых ускорений, а результаты
полученные на GPU, методами сопряженных градиентов и пополнения, эффективнее более чем в десять раз (bodyy4).
Дальнейшие развитие высопроизводительных методов обращения матриц связано с модификацией и оптимизацией алгоритма метода Шермана - Моррисона при обращении симметричных матриц, в том числе и для построения эффективных предобуславливателей на его основе для метода сопряженных градиентов.
Работа выполнена при финансовой поддержке РФФИ (грант Ц-01-31066-мол_а, 14-01-00055-а, 13-01-00101-а) и программы Президиума РАН №18 при поддержке УрО РАН (проект 12-П-1-1005).
Литература
1. Копысов, С.П., Кузьмин, И.М., Недожогин, Н.С., Новиков, А.К. Параллельные алгоритмы формирования и решения системы дополнения Шура на графических ускорителях // Ученые записки Казанского университета. Серия Физ.-мат. науки - 2012. - Т. 154. №3.
- С. 202-215.
2. Ezzatti, P., Quintana-Orti, E.S., Remon, A. Using graphics processors to accelerate the computation of the matrix inverse // J. of Supercomputing. - 2011. V.58. - P.429-437.
3. Kopysov, S.P., Kuzmin, I.M., Nedozhogin, N. S., Novikov, A.K., Sagdeeva, Y.A. Hybrid Multi-GPU solver based on Schur complement method // Lecture Notes in Computer Science - 2013.
- vol. 7979. - pp. 65-79.
4. Woodbury, M. Inverting modified matrices // Memorandum Rept. 42, Statistical Research Group, Princeton University, Princeton, NJ, 1950.
5. Sherman, J., Morrison, W.J. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix // Ann. Math. Statistics. - 1950. - V. 21. №1. - P. 124-127.
6. Фаддеев, Д.К., Фаддеева, В.Н. Вычислительные методы линейной алгебры. - СПб.: Изд-во "Лань 2002. - 736 с.
7. He, X., Holm, M., Neytcheva, M. Efficiently parallel implementation of the inverse Sherman -Morrison algorithm // Lecture Notes in Computer Science. - 2013. - V. 7782. — P. 206-219.
8. Matrix Market: URL: http://math.nist.gov/MatrixMarket/ (дата обращения: 29.08.2013)
9. Davis, T. University of Florida Sparse Matrix Collection: sparse matrices from a wide range of applications: URL: http://www.cise.ufl.edu/research/sparse/matrices/ (дата обращения:
29.08.2013)
Недожогин Никита Сергеевич, аспирант, Институт механики УрО РАН (Ижевск, Российская Федерация), [email protected].
Сармакеева Анастасия Семеновна, студент математического факультета, Удмуртский государственный университет (Ижевск, Российская Федерация), [email protected].
Копсов Сергей Петрович, д.ф.-м.н., профессор, заведующий лабораторией, Институт механики УрО РАН (Ижевск, Российская Федерация), [email protected].
Поступила в редакцию 14 марта 2014 г. 106 Вестник ЮУрГУ. Серия «Вычислительная математика и информатика»
H.C. He,n,o:*orHH, A.C. CapMaKeeBa, C.n. KonucoB
Bulletin of the South Ural State University Series “Computational Mathematics and Software Engineering”
2014, vol. 3, no. 2, pp. 101-108
SHERMAN - MORRISON HIGH-PERFORMANCE ALGORITHM FOR INVERSE MATRIX ON GPU
N.S. Nedozogin, Institute of Mechanics of Ural Branch of the RAS (Izhevsk,
Russian Federation),
A.S. Sarmakeeva, Institute of Mechanics of Ural Branch of the RAS
(Izhevsk, Russian Federation),
S.P. Kopysov, Institute of Mechanics of Ural Branch of the RAS (Izhevsk,
Russian Federation)
Matrix inversion is widely used in numerical methods, such as linear solvers, preconditioning for linear system, domain decomposition, digital image processing, etc. High-performance implementation of matrix inversion requires efficient matrix storage formats and optimal distribution of computations between computing devices. In this paper, we study the performance of traditional matrix inversion algorithms, such as LU-factorization and Gauss-Jordan, as well as the conjugate gradient method and the Sherman - Morrison formula. In the last two algorithms, matrix-vector products and scalar products are efficiently executed on multicore/manycore processors. We compare the performance of the algorithms on hybrid multi-CPU multi-GPU platforms, using the matrices from well-know test suites and from the numerical simulation of wrap spring.
Keywords: Sherman - Morrison formula, high-performance computing, inverse matrix, sparse matrix.
References
1. Kopysov S.P., Kuzmin I.M., Nedozhogin N.S., Novikov A.K. Parallelnye algoritmy formirovaniya i resheniya sistemy dopolneniya Shura na graficheskix uskoritelyax // Uchenye zapiski Kazanskogo universiteta. Seriya Fiz.-mat. nauki - 2012. - T. 154. №3. - P. 202-215.
2. Ezzatti P., Quintana-Orti E.S., Remon A. Using graphics processors to accelerate the computation of the matrix inverse // J. of Supercomputing. - 2011. V.58. - P.429-437.
3. Kopysov S.P., Kuzmin I.M., Nedozhogin N.S., Novikov A.K., Sagdeeva Y.A. Hybrid Multi-GPU solver based on Schur complement method // Lecture Notes in Computer Science -2013. - vol. 7979. - pp. 65-79.
4. Woodbury, M. Inverting modified matrices // Memorandum Rept. 42, Statistical Research Group, Princeton University, Princeton, NJ, 1950.
5. Sherman J., Morrison W.J. Adjustment of an inverse matrix corresponding to a change in one element of a given matrix // Ann. Math. Statistics. - 1950. - V. 21. №1. - P. 124-127.
6. Faddeev D.K., Faddeeva V.N. Vychislitelnye metody linejnoj algebry. - SPb.: Lan, 2002. -736 p.
7. He X., Holm M., Neytcheva M. Efficiently parallel implementation of the inverse Sherman -Morrison algorithm //Lecture Notes in Computer Science. - 2013. - V. 7782. — P. 206-219.
8. Matrix Market: URL: http://math.nist.gov/MatrixMarket/ (accessed: 29.08.2013)
9. Davis, T. University of Florida Sparse Matrix Collection: sparse matrices from a wide range of applications: URL: http://www.cise.ufl.edu/research/sparse/matrices/ (accessed:
29.08.2013)
Received 14 March 2014-
