Решение систем линейных алгебраических уравнений методом предобуславливания на графических процессорных устройствах Текст научной статьи по специальности «Математика»
Аннотация научной статьи по математике, автор научной работы — Карпенко А. П., Чернов С. К.
Рассматриваем алгоритм решения систем линейных алгебраических уравнений (СЛАУ) с предобуславливанием . Представляем параллельные алгоритмы и программы для графических процессорных устройств , реализующие основные операции этого алгоритма – операции умножения матрицы на набор векторов и решения блочно-треугольной СЛАУ. Приводим результаты широкого исследования эффективности предложенных алгоритмических и программных решений. Эти результаты показывают достаточно высокую эффективность разработанных алгоритмов и программ для умножения матрицы на набор векторов. Ускорение вычислений в этом случае составляет от четырех до шестнадцати раз. Алгоритмы и программы, предназначенные для решения блочно-треугольной СЛАУ, показали удовлетворительные результаты, которые позволяют ожидать приемлемого ускорения для практически значимых СЛАУ высокой размерности и при использовании профессиональных ГПУ.
Похожие темы научных работ по математике , автор научной работы — Карпенко А. П., Чернов С. К.
Текст научной работы на тему «Решение систем линейных алгебраических уравнений методом предобуславливания на графических процессорных устройствах»
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 • 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Решение систем линейных алгебраических уравнений
методом предобуславливания на графических процессорных
Карпенко А. П., Чернов С. К.
Россия, МГТУ им. Н.Э. Баумана apkarpenko@mail.ru lemonlime@yandex.ru
В настоящее время вычислительная гидродинамика находит широкое применение при решении, например, задач авиа-, судо - и двигателестроения. В силу высокой вычислительной сложности методов этой науки для решения указанных и других задач используют параллельные вычислительные системы различных классов и, прежде всего, кластерные системы. Для уменьшения расходов на вычисления в настоящее время в качестве альтернативы кластерам все чаще используют графические процессорные устройства (ГПУ). Статья подготовлена в рамках работ по адаптации программного комплекса FlowVision [1] российской компании «Тесис» к ГПУ
На верхнем иерархическом уровне рассматриваемая в работе задача представляет собой задачу течения жидкости или газа в некотором заданном объеме. Имеется в виду, что для описания течения используются математически модели, основанные на дифференциальных уравнениях в частных производных (ДУЧП) Навье-Стокса, Эйлера, Рейнольдса и т. д. Для численного решения этих уравнений чаще всего применяют методы конечных разностей (КР), конечных элементов (КЭ) и конечных объемов (КО).
При использовании регулярной сетки метод КР прост в реализации,
эффективен, имеет наглядную процедуру дискретизации и позволяет легко получить высокий порядок точности. Существенным недостатком метода является то, что его достоинства проявляются только при использовании регулярной сетки, которая возможна лишь в случае простой геометрии расчетной области [2]. К достоинствам метода КЭ относится простота использования его для задач с произвольной формой области решения, к недостаткам — высокая вычислительная и алгоритмическая сложность [2].
И метод КР и метод КЭ обладают еще тем недостатком, что обычно не обеспечивают в явном виде выполнения законов сохранения [3]. От этого недостатка свободен метод КО [3]. По сравнению с методов КР, метод КО обладает более высокими возможностями использования в областях произвольной формы, а по сравнению с методом КЭ - меньшими вычислительными затратами. На этом основании в программном комплексе FlowVision применяется именно метод КО.
Численное решение краевых задач для ДУЧП является весьма ресурсоемкой задачей. Так, например, современные практически значимые задачи гидрогазодинамики требуют использования неструктурированных адаптивных расчетных сеток с числом ячеек порядка 100 млн. Поэтому одним из главных направлений развития FlowVision является повышение эффективности использования вычислительных ресурсов.
В настоящее время FlowVision ориентирован на кластерные вычислительные системы, процессоры которых являются многоядерными [1]. На каждом процессоре кластера FlowVision запускает один процесс, использующий библиотеки Intel TBB и OpenMP для создания некоторого динамически определяемого числа потоков (threads). Обмен данными между этими потоками происходит через общую память процессора. Обмен данными между процессами, запущенными на разных процессорах, реализуется с помощью библиотека MPI.
В связи с тем, что в настоящее время идет активное внедрение гетерогенных кластерных вычислительных систем, в состав которых входят
ГПУ, все более актуальной становится задача адаптации программного обеспечения и, в частности, FlowVision, к таким системам. Решение этой задачи, с одной стороны, может дать заметный выигрыш в производительности, с другой стороны, сложная многоуровневая гибридная архитектура таких вычислительных систем приводит к значительному усложнению программного обеспечения.
Современными производителями ГПУ являются, прежде всего, компании AMD и NVIDIA. Основу ГПУ обоих производителей составляет набор SIMD узлов (потоковых мультипроцессоров), связанных разделяемой памятью [4]. Каждый из потоковых мультипроцессоров включает в себя несколько потоковых процессоров, блок специализированных АЛУ, обеспечивающих вычисление наиболее часто встречающихся математических функций, блок вычислений с двойной точностью, разделяемую память и т.д.
Программирование ГПУ обеих указанных компаний осуществляют с помощью универсальной технологии OpenCL [5] и технологии CUDA [6], разработанной компанией NVIDIA для ее собственных ГПУ Несмотря на универсальность технологии OpenCL, работа ориентирована на использование ГПУ NVIDIA и технологии CUDA. Данный выбор обусловлен следующими соображениями:
- неполная на момент начала работы поддержка обоими производителями стандарта OpenCL;
- более интенсивное развитие технологии CUDA по сравнению с технологией OpenCL;
- наличие у авторов опыта работы с технологией CUDA.
Компания NVIDIA предоставляет библиотеку cuBLAS, реализующую некоторые алгоритмы линейной алгебры [7]. Например, библиотека поддерживает умножение матрицы на матрицу, матрицы на вектор и т.д. Однако библиотека ориентирована только на операции с плотными матрицами и поэтому не может быть использована в нашем случае. Для работы с
разреженными матрицами предназначена ориентированная на ГПУ библиотека сыБРЛЕБЕ [8], но в этой библиотеке в настоящее время реализована операция умножения блочно-разреженной матрицы только на один вектор, в то время как в FlowVision требуется умножение такой матрицы на набор векторов.
В первом разделе работы представлена постановка задачи. Во втором разделе обсуждаем алгоритм неполных треугольных факторизаций второго порядка точности 1Ш2, распараллеливанию которого и посвящена работа. В третьем разделе рассматриваем алгоритмическую и программную реализации на ГПУ двух основных операций алгоритма 1Ш2 - операция умножения матрицы на блок векторов и операция решения блочно-треугольной системы линейных алгебраических уравнений (СЛАУ). В четвертом разделе представлены результаты исследования эффективности принятых алгоритмических и программных решений. В заключении сформулированы основные результаты работы и перспективы ее развития.
1 Постановка задачи
FlowVision позволяет производить полный комплекс действий для расчёта 3О стационарных и нестационарных течений газов и жидкостей в областях произвольной сложной формы с учетом теплообмена и химических реакций:
- импорт геометрической модели;
- автоматизированное построение расчетной сетки на основе модели;
- задание начальных и граничных условий;
- формирование математической модели в виде СЛАУ на основе метода КО;
- решение указанной СЛАУ;
- визуализация результатов моделирования.
Комплекс использует адаптивную локально измельчаемую расчетную сетку с прямоугольными ячейками и с подсеточным разрешением геометрии [9]. Последний механизм позволяет аппроксимировать
криволинейные границы без измельчения сетки и основан на том, что ячейка, через которую проходит граница, разбивается на несколько ячеек, часть из которых исключается из расчета. Оставшиеся ячейки представляют собой не прямоугольные параллелепипеды, а многогранники произвольной формы. Важно, что уравнения математической модели для этих многогранников FlowVision аппроксимирует без упрощений, что позволяет рассчитывать течения с высокой точностью даже на грубой сетке.
Метод КО, в конечном счете, сводит решение модельной краевой задачи для ДУЧП к решению СЛАУ вида
где А - (уху)-разреженная квадратная матрица, Ь - (ух 1)-вектор правой части, x - искомый (ух 1)-вектор решения. Особенностями СЛАУ (1)
являются высокая размерность и возможная плохая обусловленность.
Для решения СЛАУ разработано большое число различных методов, наиболее простыми из которых являются прямые методы (метод исключения Гаусса, Гаусса-Жордана и прочие). Хорошо известно, что методы этого класса не эффективны при решении плохо обусловленных систем и систем высокой размерности. Как правило, в вычислительной практике применяют итерационные методы и, в частности, методы, основанные на предобуславливателях [10].
Известно значительное число методов предобуславливания: приближенные обратные методы; приближенные факторизованные обратные методы; методы, основанные на неполной треугольной факторизации и т. д. Методы различаются затратами на построение предобуславливателя, свойствами сходимости и масштабируемости при параллельной реализации.
В комплексе FlowVision используется алгоритм неполных треугольных факторизаций второго порядка точности ^Ш, являющийся несимметричным обобщением алгоритма ЮШ [10]. Алгоритм показал высокую
вычислительную эффективность для широкого класса задач, в том числе, для задач вычислительной гидрогазодинамики [11 ]. Ставится задача реализации этого алгоритма на ГПУ.
2. Распараллеливание алгоритма 1ЬШ
Процесс факторизации в соответствии с алгоритмом КШ основан на матричном соотношении
A + E = LU + LR + WU, где A - масштабированная матрица системы (1), L и U — нижняя и верхняя треугольные части элементов первого порядка точности, W и R -аналогичные части элементов второго порядка точности, Е - матрица ошибок.
В итерационной схеме алгоритма итерации производятся с
матрицей . Алгоритм требует хранения в памяти ЭВМ всех
элементов матриц А, L, и, небольшой части элементов матриц W и R и включает в себя следующие операции [11]:
1) построение предобуславливателя (неполная факторизация матрицы),
2) умножение матрицы на вектор,
3) решение систем уравнений с верхней и нижней треугольными матрицами,
4) векторные операции.
Для обеспечения высокой эффективности реализации на ГПУ алгоритма ^и2 необходимо согласованное распараллеливание всех указанных операций. Рассмотрим первые три из них.
Построение предобуславливателя. В блочном виде рекуррентное выражение для операции треугольной факторизации имеет вид
Аи а1,2 А1,3 L1,1 0 0 иц и1,2 и1,3
а2,1 а2,2 а2,3 = 12,1 12,2 0 X 0 и 2,2 и 2,3
А3,1 а3,2 А3,3 _ ^,1 13,2 ^'3,3 0 0 и3,3
Здесь блоки а22, 12 2, и2 2 имеет единичный размер; блоки а21, а23, /21, и23 являются строками; блоки а12, а32, /32, и12 - столбцами; Li , и< ^ - квадратные матрицы. На каждой итерации необходимы следующие операции:
1) на основе уравнения l22 u22 = a22 -l21 u12 вычислить значения
компонентов строк l2 2, u
2) из уравнений u
компоненты строки и23 и столбца 13 2 .
Зависимости по данным для строки матрицы и представленного алгоритма иллюстрирует рисунок 1. Для столбца матрицы Ь зависимость по данным получается симметричным отражением представленных зависимостей относительно главной диагонали матрицы.
Рисунок 1 - Зависимости по данным при вычислении компонентов строки
матрицы и в процессе факторизации
Различные варианты метода неполных разложений, в том числе алгоритмы 1СН2, 1Ьи2, представляют собой некоторые модификации шагов 1, 2 [11].
Умножение матрицы на вектор осуществляется по формуле
Z A .x. , i = 0,1. (v — 1):
из которой следует, что при выполнении данной операции рекуррентная зависимость по данным отсутствует.
Решение систем уравнений. Рекуррентные уравнения для решения треугольных систем с матрицами L и U имеют, соответственно, вид