Исследование эффективности переупорядоченного метода BiCGStab на вычислительных системах СКИФ МГУ «Чебышёв» и «Ломоносов» Текст научной статьи по специальности «Математика»

УДК 519.852.67

ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ ПЕРЕУПОРЯДОЧЕННОГО МЕТОДА BICGSTAB НА ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ СКИФ МГУ «ЧЕБЫШЁВ» И «ЛОМОНОСОВ»

Б.И. Краснопольский

AN INVESTIGATION OF THE EFFICIENCY OF REORDERED BICGSTAB METHOD ON COMPUTER SYSTEMS SKIF MSU «CHEBYSHEV» AND «LOMONOSOV»

B.I. Krasnopolsky

В представленной работе обсуждаются результаты исследования эффективности и масштабируемости переупорядоченного итерационного метода BiCGStab на вычислительных системах СКИФ МГУ «Чебышёв» и «Ломоносов». Основное внимание уделяется решению больших систем линейных алгебраических уравнений с сильно-разреженной матрицей, возникающих при решении уравнения Пуассона в различных геометрических областях. Рассматривается вопрос целесообразности применения гибридных моделей программирования для решения задач линейной алгебры на многопроцессорных вычислительных системах с распределенной памятью.

Ключевые слова: итерационные методы, BiCGStab, разреженные матрицы, масштабируемость, уравнение Пуассона.

The results of the efficiency investigation of reordered BiCGStab method on computer systems SKIF MSU «Chebyshev» and «Lomonosov» are presented. The main attention is focused on solving the systems of linear algebraic equations with large sparse matrices for the Poisson equation in different geometries. The reasonability of usage the hybrid programming models for the linear algebra problems on distributed memory computer systems is discussed.

Keywords: iterative methods, BiCGStab, sparse matrices, scalability, the Poisson equation.

Введение

Одним из наиболее трудоемких этапов численного решения системы уравнений гидродинамики сеточными методами является решение уравнения неразрывности, сводящегося к уравнению Пуассона для давления. Зачастую этот этап оказывается доминирующим и занимает подавляющую часть времени (более 90%). Для простых геометрий расчетных областей и ортогональных систем координат, а так же областей, отображаемых преобразованиями систем координат в прямоугольные области, разработаны прямые быстрые методы решения уравнения Пуассона (например, [1]). Однако, для сколь-нибудь сложных конфигураций расчетных областей, а так же неортогональных расчетных сеток данные методы оказываются

неприменимыми. В таком случае достаточно популярным является подход, при котором используется итерационный метод градиентного типа, например метод ВЮС81аЬ [2], с предо-буславливателем для ускорения сходимости. В качестве предобуславливателей долгое время использовалось частичное Ы1-разложение [3] с различной степенью заполненности. Однако, при всех своих плюсах (в первую очередь, простота и эффективность при последовательной реализации), данный метод обладает рядом значительных недостатков. Основные из них -это ухудшение эффективности предобуславливания по мере увеличения матрицы системы и высокая чувствительность к переупорядочиванию элементов матрицы. Последний аспект является критическим при параллельных вычислениях, поскольку эффективное распределение матрицы системы между вычислительными процессами также требует переупорядо-чивания элементов матрицы. Ввиду этого для решения уравнений эллиптического типа в произвольных областях набирают популярность многосеточные методы [4], которые могут использоваться как отдельно взятые решатели, так и эффективные предобуславливатели.

1. Переупорядоченный метод BiCGStab

Алгоритм классического метода BiCGStab [2] приведен на рис. 1. Каждая итерация этого метода включает в себя две операции умножения разреженной матрицы на вектор (и две операции предобуславливания для предобусловленного варианта метода), четыре скалярных произведения векторов и двенадцать операций сложения векторов/умножения на число. При традиционном подходе к распараллеливанию итерационных градиентных методов с сильно-разреженной матрицей, исходная матрица системы и вектора распределяются между вычислительными процессами построчно [5], и каждый из локальных фрагментов матрицы также разбивается на блоки пропорционально количеству строк, пришедшихся на каждый из вычислительных процессов. Такая схема распределения данных между вычислительными процессами позволяет эффективно организовать вычисление операций сложения векторов/умножения на число и умножение матрицы на вектор: операции с векторами выполняются локально, тогда как время на пересылку данных для вычисления произведения матрицы на вектор может быть эффективно «замаскировано» вычислениями с локальным фрагментом вектора. Наибольшие сложности возникают при распараллеливании скалярных произведений векторов, поскольку после вычисления скалярного произведения локального фрагмента вектора требуется взаимодействие всех процессов для обмена данными каждого вычислительного процесса с каждым. Такое взаимодействие может быть реализовано как одна коллективная операция или набор коммуникаций «точка-точка» для каждого процесса с каждым. Второй вариант может оказаться предпочтительным, если алгоритмически имеется возможность организовать порядок вычислений так, что время на выполнение этих коммуникаций будет «замаскировано» другими вычислениями. В противном случае операция скалярного произведения векторов является точкой глобальной синхронизации всех вычислительных процессов, что нежелательно при реализации параллельной программы.

Приведенный на рис. 1 классический алгоритм метода BiCGStab демонстрирует строгую последовательность вычислений, то есть для выполнения вычислений г-й строки алгоритма необходимо, чтобы были завершены вычисления всех предыдущих г — 1 строк. Таким образом, присутствие нескольких скалярных произведений векторов, которые не могут быть алгоритмически переставлены, а время на коммуникации для завершения этих операций «замаскировано» другими вычислениями, приводит к наличию трех точек глобальной синхронизации вычислений, и, как следствие, к ухудшению эффективности параллельной реализации исходного метода.

Известно несколько работ [6, 7], посвященных попыткам переформулировать непредо-бусловленный метод BiCGStab с целью уменьшить количество операций скалярного про-

1. х0 = начальное приближение; г0 = Ь — Аха

2. ро = (го, г0)

3. ро =

для ; = 0,1,...

4. и = Ар,,-

5- а> =

6. зу = rj — ар

7. t =

О , . _ >3] )

шз (г,£)

9. х]+1 = + азР] + шг5^; проверка сходимости

10. Гу+1 = Эу — со^

11. рт = (г;+1,г0)

12. Я. = £Ш.2£

гз

13. р^+1 = г,-+1 + ^{р5 - шр)

конец Рис. 1. Классический метод В1ССЭ1аЬ

изведения векторов, приводящих к глобальной синхронизации вычислительных процессов. В [6] предложен вариант метода с двумя глобальными точками синхронизации при вычислении скалярных произведений векторов, а в [7] позднее был предложен вариант с одной точкой синхронизации. Основываясь на идеях и подходах, изложенных в [6, 7], в [8, 9] для предобусловленного варианта метода, была предложена модификация, которая позволяет избежать глобальной синхронизации всех вычислительных процессов на каждой итерации. Предложенный метод имеет только четыре дополнительных операции сложения векторов/умножения на число по сравнению с исходным алгоритмом, в то время как IBiCGStab [7] требует два дополнительных скалярных произведения и шесть операций с векторами. Скорость сходимости модифицированного метода идентична скорости сходимости предобусловленного классического метода, поскольку в предложенной модификации изменена только последовательность вычисления некоторых векторов. Отличительной же чертой этого метода является полное отсутствие точек глобальной синхронизации процессов. Скалярные произведения векторов сгруппированы в три блока, причем для каждого из них имеется возможность «замаскировать» время обмена данными вычислением предобуславливания, которое зачастую является наиболее трудоемкой операцией на итерации.

2. Особенности параллельной реализации

Исследование эффективности предложенной в [8, 9] модификации алгоритма было проведено на вычислительных системах СКИФ МГУ «Чебышёв» и «Ломоносов». Для этого была разработана гибридная параллельная программа, реализующая предложенный метод с алгебраическим многосеточным предобуславливателем (для построения предобуславливания использована свободно-распространяемая библиотека численных методов hypre [10]). Взаимодействие между узлами было реализовано через MPI, тогда как обмен данными между процессами внутри одного узла реализован через POSIX Shared Memory. Для хранения исходной и промежуточных матриц использованы специализированные форматы CRS и DMSR (например, [11]).

Обычно, для построения оптимального переупорядочивания применяются методы раз-

биения графов или гиперграфов, которые строятся по исходной матрице системы, подлежащей переупорядочиванию. Был проведен ряд исследований по изучению эффективности различных алгоритмов, реализованных в свободно-распространяемых библиотеках METIS, hMETIS, ParMETIS, Scotch, PaToH и ряде других. По результатам этих исследований были выбраны библиотека Scotch и Pt-Scotch для разбиения локального и распределенного графов соответственно.

3. Тестовые задачи

Тестовыми задачами для исследуемого метода и его реализации были выбраны:

1. задача Неймана для уравнения Пуассона с постоянной правой частью в кубической области на равномерной сетке;

2. задача Неймана для уравнения Пуассона для давления с постоянной правой частью в расчетной области, представляющей собой канал с препятствием на стенке в виде куба (задача обтекания куба турбулентным потоком [12]; 5,45 млн. ячеек).

При проведении расчетов критерием завершения итерационного процесса было задано условие1:

llsj+i ~ xi\\ < 10-б_

1М

В качестве характеристики эффективности параллельной программы использовано ускорение решения задачи, определяемое как отношение времени решения задачи на одном вычислительном ядре (процессоре) ко времени решения на заданном количестве вычислительных ядер (процессоров):

Тх

Pn = tn

4. Результаты

Первоначально была реализована параллельная программа с взаимодействием вычислительных процессов только через MPI. На рис. 2 представлены полученные графики зависимости ускорения решения задачи от количества используемых вычислительных узлов для уравнения Пуассона в кубической расчетной области размером 8 млн. ячеек (количество ядер в легенде графика соответствует количеству вычислительных процессов, запущенных на каждом из узлов вычислительной системы). Следует отметить относительно слабые результаты, полученные для большого числа вычислительных процессов на узлах. Выигрыш от использования всех 8 вычислительных ядер на каждом узле в сравнении с 4 ядрами, при количестве используемых узлов до 16, оказывается крайне незначительным, тогда как при большем количестве узлов использование 8 ядер на узле оказывается менее эффективным. Максимальное ускорение в такой конфигурации не превышает 35 раз. Эффективность использования 4 ядер с каждого узла оказывается несколько выше: насыщение достигается на отметке 48 узлов (196 вычислительных ядер) при ускорении в 66 раз. Наиболее удачной оказывается конфигурация, при которой задействованы лишь два вычислительных ядра с каждого узла: на малом количестве узлов преимущество от использования большего количества ядер оказывается незначительным, а насыщение коммуникационной сети, соединяющей вычислительные узлы, наступает значительно позже. В такой конфигурации

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

160................

♦ 1 ядро/узел

количество узлов

Рис. 2. Ускорение решения системы уравнений в зависимости от количества узлов и вычислительных процессов на каждом из узлов (Ьурге 2.4.0Ь, СКИФ МГУ «Че-бышёв»)

180...... .......................................... ■ ■ ■ .......

1 тред

количество узлов

Рис. 3. Ускорение решения системы уравнений в зависимости от количества узлов и вычислительных процессов, объединенных через общую память на каждом из узлов (hypre 2.4.0b, СКИФ МГУ «Чебышёв»)

получено ускорение до 140 раз на 112 вычислительных узлах (224 ядра), причем до этой отметки зависимость ускорения от количества узлов носит практически линейный характер. Для каждого из вариантов оптимальным оказывается использование порядка 200 - 250 вычислительных процессов, что эквивалентно локальной области около 30 - 35 тыс. строк исходной матрицы, приходящейся на один вычислительный процесс.

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

количество узлов

Рис. 4. Ускорение решения системы уравнений в зависимости от количества узлов и размера матрицы относительно времени расчета на 4 узлах (32 вычислительных ядра, Ьурге 2.6.0Ь, СКИФ МГУ «Чебышев»)

улучшить ситуацию с использованием всех имеющихся ядер в вычислительных узлах и уменьшить нагрузку на коммуникационную сеть.

Такая версия программы была исследована на задаче решения уравнения Пуассона в канале с кубическим препятствием [12]. Соответствующие графики зависимости ускорения от количества вычислительных узлов приведены на рис. 3 (количество тредов в легенде графика соответствует количеству МР1-процессов внутри узла, взаимодействие между которыми реализовано через общую память). Сравнивая рис. 2 и 3, можно отметить, что для второй тестовой задачи выигрыш от использования всех имеющихся восьми ядер с каждого узла при использовании только МР1-обменов (кривая на графике для одного треда) оказывается еще меньшим, чем в предыдущем случае, и ускорение не превышает 20 раз. Этот результат объясняется уменьшением локальной области, приходящейся на каждый вычислительный процесс. Размер локального фрагмента матрицы при пиковом значении ускорения уменьшается с 35 до 28 тыс. строк.

Вместе с тем, представленный график наглядно демонстрирует преимущество гибридного подхода к реализации алгоритма, который позволяет существенно снизить нагрузку на коммуникационную сеть кластера и добиться эффективной утилизации всех имеющихся процессорных ядер на узлах. Пиковое ускорение для данной тестовой задачи составило 156 раз (0,77 секунды) в сравнении с последовательно исполняемой программой (120,1 секунды) при использовании 224 узлов (1792 ядра). Локальная область данных исходной матрицы, приходящаяся на каждый процесс, составила примерно 3 тыс. строк, что является крайне малым показателем для таких задач.

Далее была исследована зависимость масштабируемости решателя от размера матрицы системы уравнений. Такие расчеты были проведены для матриц размером 8, 27, 42,9 и 64 млн. строк2. Соответствующие графики зависимости ускорения от количества вычислительных узлов приведены на рис. 4. Поскольку матрицы размером более 8 млн. строк

23дееь и далее в расчетах использована обновленная версия библиотеки Ьурге 2.6.0Ь, в которой был расширен набор параметров для построения алгебраического многосеточного предобуславлива-теля.

Таблица 1

Пиковые значения ускорения для различных тестовых матриц (Ьурге 2.6.0Ь, СКИФ МГУ «Чебышёв»)

Размер матрицы, млн.строк Оптимальное количество узлов Пиковое ускорение

8 192 195

27 192 336

42,9 192 340

64 192 375

не умещаются в оперативную память одного вычислительного узла, то приведенные ускорения нормированы на время решения задачи на 32 ядрах (4 вычислительных узла). Для тестовых матриц 27 млн. строк и более приведенные графики ускорения носят линейный характер вплоть до 192 узлов. При увеличении размера тестовой матрицы наблюдается некоторое улучшение показателей масштабируемости решателя, что обусловлено увеличением размеров локальных фрагментов данных каждого из процессов. Поведение кривых при увеличении количества использованных узлов до 224 несколько странно, и, вероятно, вызвано особенностями работы аппаратной платформы во время проведения тестов.

Основываясь на практически линейном поведении кривых ускорения на малом количестве узлов и приняв за исходные результаты масштабируемости для матрицы размером 8 млн. строк, можно получить оценочные значения ускорения относительно последовательно исполняемых вариантов. В этом случае коэффициент ускорения определяется как

О T‘i2 Tf Т32

■nv = TfTTfT =

i32 iN J-N

где коэффициент К32 для матрицы 8 млн. строк оказался равен Ку2 = 9,27. Полученные оценочные значения ускорения приведены в табл. 1. С учетом того, что при увеличении размеров матрицы, наблюдается улучшение показателей масштабируемости на фиксированном количестве узлов, можно утверждать, что такой подход позволяет получить оценку снизу реального уровня масштабируемости решателя.

Сравнительные графики ускорения решения системы уравнений для матрицы размером 8 млн. строк на вычислительных системах СКИФ МГУ «Чебышёв» и «Ломоносов» приведены на рис. 5. Данные, полученные на вычислителной системе «Ломоносов», демонстрируют существенно лучшие показатели масштабируемости на малом количестве узлов, ввиду, в первую очередь, более эффективного доступа к памяти и работы ядер внутри вычислительного узла. Ускорение в пределах одного узла, в среднем, составляет 3,2 раза в сравнении с 2,4 для «Чебышёва». Вместе с тем, пиковые значения ускорения для обеих вычислительных систем оказываются примерно одинаковыми. Незначительное преимущество «Чебышёва» может быть вызвано различиями в заданных параметрах предобуславлива-теля и использованием различных реализаций библиотеки MPI. Однако, если принять во внимание тот факт, что время решения системы уравнений на одном вычислительном ядре системы «Ломоносов» оказывается меньшим в 1,7 раза, то наблюдаемое насыщение кривых на меньшем количестве узлов представляется закономерным: время проведения вычислений с локальными фрагментами данных, сосредоточенными в каждом вычислительном процессе, для «Ломоносова» оказывается существенно меньшим. С этой точки зрения, для сравнения вычислительных систем более показателен третий из приведенных на рис. 5 графиков. Этот график отражает ускорение решения задачи относительно времени решения на одном ядре «Чебышёва» и демонстрирует более чем 300-кратное ускорение.

количество узлов

Рис. 5. Ускорение решения системы уравнений в зависимости от количества узлов на различных аппаратных платформах (матрица 8 млн. строк, Ьурге 2.6.0Ь)

Таблица 2

Пиковые значения ускорения для различных тестовых матриц (Ьурге 2.6.0Ь, «Ломоносов»)

Размер матрицы, млн. строк Оптимальное количество узлов Пиковое ускорение

8 96 178

27 224 337

64 320 522

125 416 738

216 608 1134

343 768 1417

Полученные на вычислительной системе «Ломоносов» графики зависимости масштабируемости решателя от размера тестовой матрицы приведены на рис. 6. Как и для «Че-бышёва», ускорение решения системы уравнений было определено относительно времени решения на 32 узлах (256 вычислительных ядер). Приведенный график демонстрирует близкую к линейной масштабируемость на больших матрицах вплоть до 512 узлов (4096 вычислительных ядер) и положительный эффект от увеличения задействованных ресурсов вплоть до 768 узлов (6144 вычислительных ядра). Аналогичным образом получены оценочные значения абсолютных показателей масштабируемости метода и реализованного решателя, которые приведены в табл. 2, где соответствующий коэффициент -¿^256 Для матрицы 8 млн. строк составил К^56 — 95, 93. Результаты пикового ускорения показывают практически линейный рост относительно размера матрицы системы.

Времена решения различных систем уравнений и соответствующее количество итераций до сходимости итерационного процесса приведены в табл. 3. Реализованный решатель демонстрирует слабую зависимость скорости сходимости как от размера матрицы, так и от количества вычислительных ядер, что обеспечивает хорошие результаты масштабируемости при решении различных систем линейных алгебраических уравнений.

количество узлов

Рис. 6. Ускорение решения системы уравнений в зависимости от количества узлов и размера матрицы относительно времени расчета на 32 узлах (256 вычислительных ядер; Ьурге 2.6.0Ь, «Ломоносов»)

Таблица 3

Время решения системы линейных алгебраических уравнений для различных тестовых матриц при использовании 128 узлов (1024 вычислительных ядра; Ьурге 2.6.0Ь, «Ломоносов»)

Размер матрицы, млн. строк Количество итераций Время решения, сек

8 10 0,19

27 11 0,4

64 12 0,88

125 13 1,69

216 14 3,13

343 14 5

Заключение

В работе представлены результаты исследования масштабируемости предобусловлен-ного переупорядоченного метода BiCGStab для решения больших систем линейных алгебраических уравнений с сильно разреженной матрицей. Благодаря гибридному подходу (МР1+8ЬМ), реализованный метод с алгебраическим многосеточным предобуславливанием демонстрирует хорошие результаты масштабируемости (более 3 порядков) для тестовых задач на большом количестве вычислительных узлов: до 200 узлов на вычислительной системе СКИФ МГУ «Чебышёв» и до 1000 узлов на вычислительной установке «Ломоносов». Реализованный решатель демонстрирует слабую зависимость скорости сходимости итерационного процесса как от размера матрицы, так и от количества вычислительных ядер, что обеспечивает хорошие результаты масштабируемости при решении различных систем линейных алгебраических уравнений.

Представленная работа частично поддержана грантом РФФИ №09-08-00390-а.

Приведенные в работе результаты исследований были получены на вычислительных системах СКИФ МГУ « Чебышев» и «Ломоносов> суперкомпъютерного комплекса МГУ им. М.В. Ломоносова. Автор благодарит Вл.В. Воеводина и С.А. Жуматия за предоставленную возможность проведения тестов на указанных вычислительных ресурсах и оперативную помощь при решении технических вопросов.

Статья рекомендована к печати программным комитетом Международной суперком-пьютерной конференции <Научный сервис в сети Интернет: суперкомпьютерные центры и задачи».

Литература

1. Swarztrauber, P.N. A direct method for the discrete solution of separable elliptic equations / P.N. Swarztrauber // SIAM J. on Numerical Analysis. - 1974. - Vol. 11. -P. 1136 - 1150.

2. Van der Vorst, H.A. BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems / H.A. van der Vorst // SIAM J. on Scientific and Statistical Computing. - 1992. - Vol. 13. - P. 631 - 644.

3. Saad, Y. Iterative methods for sparse linear systems, 2nd edition / Y. Saad. - SIAM, 2003.

4. Trottenberg, U. Multigrid / U. Trottenberg, C.W. Oosterlee, A. Schuller. - N. Y.: Academic Press, 2001.

5. Arbenz, P. Introduction to Parallel Computing - A practical guide with examples in C. / P. Arbenz, W. Petersen. - Oxford Texts in Applied and Engineering Mathematics, №. 9. Oxford University Press, 2004.

6. Jacques, T. Electromagnetic scattering with the boundary integral method on MIMD systems / T. Jacques, L. Nicolas, C. Vollaire // High-Performance Computing and Networking, Lecture Notes in Computer Science. - 1999. - Vol. 1593. - P. 1025 - 1031.

7. Yang, L. The improved BiCGStab method for large and sparse unsymmetric linear systems on parallel distributed memory architectures / L. Yang, R. Brent // Proceedings of Fifth International Conference on Algorithms and Architectures for Parallel Processing (ICA3PP’02). - 2002. - P. 324 - 328.

8. Краснопольский, В.И. Модифицированный метод BiCGStab для высокопроизводительных вычислений / Б.И. Краснопольский // Конференция молодых ученых Института механики МГУ им. Ломоносова (в печати).

9. Krasnopolsky, В. The reordered BiCGStab method for distributed memory computer systems / B. Krasnopolsky // Procedia Computer Science. - 2010. - Vol. 1. - P. 213 -218.

10. Hypre: a library of high performance preconditioners /

https: / / comput ation. llnl .gov/ease/linear _ solvers / sis _ hypre. html

11. Tuminaro, R.S. Official Aztec Users Guide, Version 2.1 / R.S. Tuminaro, M.A. Heroux, S.A. Hutchinson, J.N. Shadid. - 1999.

12. Direct numerical simulation of turbulent flow around a wall-mounted cube: spatio-temporal evolution of large-scale vortices / A. Yakhot, T. Anor, H. Liu, N. Nikitin // Journal of Fluid Mechanics. - 2006. - Vol. 566. - P. 1 - 9.

Борис Иосифович Краснопольский, кандидат физико-математических наук, научный сотрудник лаборатории общей аэродинамики НИИ механики МГУ, krasnopolsky @imec. msu. ru.

Поступила в редакцию 29 ноября 2010 г.