Алгоритмические особенности создания многосеточного решателя СЛАУ на вычислительных системах с графическими ускорителями Текст научной статьи по специальности «Математика»

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

УДК 519.612.2

АЛГОРИТМИЧЕСКИЕ ОСОБЕННОСТИ СОЗДАНИЯ МНОГОСЕТОЧНОГО РЕШАТЕЛЯ СЛАУ НА ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ С ГРАФИЧЕСКИМИ УСКОРИТЕЛЯМИ

© 2014 г. Б.И. Краснопольский,1 А.В. Медведев2

1НИИ механики МГУ им. М.В. Ломоносова 2ЗАО «Т-Сервисы», Москва

krasnopolsky@imec.msu.ru

Поступила в редакцию 07.11.2013

На практических примерах проводится сопоставление многосеточных методов SAMG и САМО для решения систем линейных алгебраических уравнений, полученных при дискретизации эллиптических уравнений. Обсуждаются ключевые особенности и проблемы реализации итерационных методов подпространства Крылова и многосеточных методов на системах с графическими ускорителями (ти1ШРЦ).

Ключевые слова: эллиптические уравнения, системы линейных алгебраических уравнений, многосеточные методы, графические ускорители.

Введение

Решение систем линейных алгебраических уравнений долгие годы является одной из типичных задач для вычислительных систем. Такая востребованность обусловлена широким спектром проблем, сводящихся в конечном счете к решению систем уравнений: например, большинство задач механики сплошных сред при использовании разностных методов и аппроксимации уравнений приводит к большим сильно-разреженным системам линейных алгебраических уравнений. В зависимости от типа исходного дифференциального уравнения и способа его аппроксимации могут использоваться различные численные методы для решения таких систем. К методам наиболее общего вида, активно используемым в настоящее время, можно отнести итерационные методы подпространства Крылова [1] в сочетании с методами неполной факторизации [1] или многосеточными методами [2] в качестве предобуслав-ливателей. При этом многосеточные методы, помимо хорошего потенциала масштабируемости, оказываются предпочтительнее ввиду ряда ключевых особенностей, как то: экономичная асимптотика общего количества операций 0(Щ, где N - размер матрицы системы; независимость скорости сходимости от размера задачи; инвариантность относительно переупорядочения эле-

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

Одной из целей проводимой работы является оценка перспектив применения графических ускорителей (GPU) для решения больших разреженных систем линейных алгебраических уравнений, и в частности, для пакета OpenFOAM применительно к задачам судостроения. Актуальность данного вопроса обусловлена потребностью проведения длительных расчетов нестационарных течений в областях сложных геометрий на подробных сетках (100 млн. ячеек и более). Отсутствие гибридных реализаций методов в OpenFOAM (MPI+OpenMP, MPI+Pthreads и пр.) существенным образом ухудшает потенциал масштабируемости приложения при использовании всех ресурсов современных многоядерных процессоров. Одним из альтернативных вариантов ускорения расчетов может явиться использование различных типов ускорителей (GPU или Intel Xeon Phi). Подробный обзор проектов, связанных с реализацией многосеточных методов на графических ускорителях, был представлен в [3]. Отметим лишь, что ни одна из доступных на данный момент реализаций многосеточных методов не обеспечивает стабильной возможности эффективного использо-

Рис. 1. Топология расчетных сеток тестовых задач моделирования течения в ступени погружного насоса

Таблица 1

Характеристики аппаратных платформ_______________________________

Название «Ломоносов» «Нижеголь» «Зилант»

Тип процессора Intel Nehalem-EP Intel Sandy Bridge-EP AMD Magny-Cours

Модель процессора X5570 (2x4 ядра) E5-2665 (2x8 ядер) Opteron 6174 (2x12 ядер)

Кол-во карт GPU/узел 2 2 -

Модель GPU NVIDIA Fermi X2070 NVIDIA Fermi M2090 -

Интерконнект IB QDR IB QDR IB DDR

вания распределенных вычислительных систем с графическими ускорителями (multiGPU).

Многосеточные методы

Алгебраические и геометрические многосеточные методы являются тематикой активных исследований на протяжении нескольких десятков лет. За это время разработано большое количество модификаций методов, ориентированных на различные области применения: решение систем линейных алгебраических уравнений для эллиптических дифференциальных уравнений с гладкими и негладкими функциями, решение систем дифференциальных уравнений, решение задач пластических деформаций материалов и другие. Краткий обзор таких методов приведен, например, в [4].

Для решения эллиптических уравнений из всего многообразия алгебраических методов обычно используются два: классический алгебраический многосеточный метод (classical algebraic multigrid, CAMG) [2] или многосеточный метод, основанный на агрегировании элементов (smoothed aggregation multigrid, SAMG) [5]. В зависимости от вида матрицы, а также необходимости перестраивать иерархию матриц при каждом следующем решении системы уравнений, преимущественным может оказаться использование как одного, так и другого метода. Темой данной статьи является сравнение двух типов многосеточных методов, CAMG и SAMG, на репрезентативных тестовых задачах, а также оценка перспектив использования этих методов при проведении расчетов на графических ускорителях.

Сравнение многосеточных методов

Тестовые заВачи. Результаты сравнения рассматриваемых численных методов обсуждаются на примере решения двух систем уравнений. Данные системы уравнений получены из пакета OpenFOAM в ходе расчета течения в ступени погружного насоса; матрицы соответствуют разностной аппроксимации эллиптического уравнения для давления в рабочем колесе и направляющем аппарате. Размеры многоблочных структурированных сеток для указанных деталей составили 2.7 млн. и 4.9 млн. неизвестных соответственно. Общий вид расчетных областей и топология сеток представлены на рис. 1.

Для решения тестовых матриц в работе использовался итерационный метод подпространства Крылова BiCGStab [6] с многосеточным методом в качестве предобуславливателя. В ходе предобуславливания выполнялась одна итерация V-цикла многосеточного метода. В случае классического многосеточного метода для построения иерархии матриц была использована библиотека hypre [7]. Построение иерархии матриц для SAMG было выполнено по аналогии с методами, реализованными в библиотеке CUSP [8, 9]. В обоих случаях использовался набор параметров методов, обеспечивающих минимальные времена решения имеющихся систем уравнений. На этапе решения в качестве пре- и пост- сглаживателя применялся итерационный метод Чебышева. Для CAMG выполнялись две итерации, для SAMG приведено сравнение результатов для двух и трех итераций.

Характеристики аппаратных платформ, использованных в ходе тестов, приведены в табл. 1.

Таблица 2

Иерархия матриц многосеточных методов, матрица 2.7 млн. неизвестных

Уровень Размер матрицы Количество ненулевых элементов Среднее количество элементов в строке

CAMG SAMG CAMG SAMG CAMG SAMG

0 2712320 2712320 18663280 18663280 6.88 6.88

1 1078627 364846 16215587 12471418 15.03 34.18

2 427949 20661 12282755 2250385 28.70 108.92

3 181283 1132 8280097 114926 45.67 101.52

4 75393 107 4617527 7487 61.25 69.97

5 28974 24 2079474 576 71.77 24.00

6 10320 - 757444 - 73.40 -

7 3464 - 230534 - 66.55 -

8 1123 - 60759 - 54.10 -

9 356 - 13578 - 38.14 -

10 109 - 2675 - 24.54 -

11 46 - 798 - 17.35 -

Иерархия матриц. Основные параметры иерархии матриц многосеточных методов для тестовой задачи размером 2.7 млн. неизвестных приведены в табл. 2; для второй тестовой задачи общие тенденции и закономерности оказываются аналогичными. Одним из главных различий между двумя рассматриваемыми методами является существенная разница в количестве уровней порождаемой иерархии матриц. CAMG несмотря на использование агрессивного подхода при выборе элементов матрицы следующего уровня (aggressive coarsening) демонстрирует не слишком высокую скорость агрегирования элементов: при переходе между уровнями размер матриц изменяется в среднем в 2.5-3 раза. Изменением значений ряда параметров метода возможно добиться более высокой скорости агрегирования элементов, однако следствием этого оказывается существенное ухудшение качества предобуславливателя. Это приводит к общему замедлению скорости сходимости метода и увеличению суммарного времени решения задачи. Для второго метода, SAMG, указанный показатель оказывается более предпочтительным, размер матриц при переходе между уровнями меняется почти на порядок. Благодаря этому иерархия матриц суммарно содержит также вдвое меньшее количество ненулевых элементов, что обеспечивает заметно меньший объем вычислений для матричных операций и

позволяет сократить затраты оперативной памяти на хранение матриц.

Скорость сходимости, времена решения и масштабируемость решателя. Использованный в расчетах набор параметров методов выбирался таким образом, чтобы обеспечить близкие к оптимальным времена решения тестовых систем уравнений. При этом анализ скорости сходимости относительной невязки показал, что рассматриваемые методы обеспечивают близкие скорости сходимости, а разброс для достижения какого-либо порогового значения не превышает нескольких итераций. Как и можно было предполагать, использование одной дополнительной итерации метода Чебышева в качестве сглаживателя обеспечивает несколько лучшую скорость сходимости. Полученные графики зависимости скорости сходимости от количества итераций методов приведены на рис. 2. В дополнение на графиках отмечены два характерных для CFD-расчетов, в том числе при решении СЛАУ для давления, уровня значений относительной невязки 10-6 и 10-8. Времена решения СЛАУ и количество итераций для трех исследуемых комбинаций методов приведены в табл. 3. Так, время расчета одной итерации метода BiCGStab с классическим алгебраическим многосеточным предобуславливателем оказывается в полтора раза больше по сравнению с SAMG(2), хотя количество операций сглажива-

Номер итерации Номер итерации

Рис. 2. Графики сходимости относительной нормы невязки. Слева - матрица 2.7 млн. неизвестных; справа -матрица 4.9 млн. неизвестных

Таблица 3

______Времена решения, с, и количество итераций до сходимости для тестовых СЛАУ__________________

Матрица Рабочее колесо, 2.7 млн. неизвестных Направляющий аппарат, 4.9 млн. неизвестных

Метод САМО SAMG(2) SAMG(3) САМО SAMG(2) SAMG(3)

Время 1 итерации 3.45 2.23 2.88 6.16 4.02 5.23

Кол-во итераций (10-6) 11 13 10 12 17 14

Время решения (10-6) 38.0 29.0 28.8 73.9 68.3 73.2

Кол-во итераций (10-8) 14 17 13 16 19 17

Время решения (10-8) 48.3 37.9 37.4 98.6 76.4 88.9

ния и прочие параметры решателя в данном случае идентичны. Такая разница обусловлена вдвое большим количеством уровней и суммарным количеством ненулевых элементов иерархии матриц, порождаемой методом САМС (табл. 2). Увеличение количества итераций сглаживания на каждом уровне многосеточного метода с двух до трех приводит к росту времени расчета на 30%, хотя таковое все равно на 20% меньше времени расчета одной итерации САМС.

Наилучшие результаты по скорости сходимости до достижения двух пороговых значений демонстрируют предобуславливатели САМС и SAMG(3), при этом разброс не превышает двух итераций в пользу того или иного метода. Более «легковесный» метод SAMG(2) демонстрирует худшую скорость сходимости и затрачивает на 3-5 итераций больше по сравнению с оптимальным методом. Однако если сравнить время решения тестовых СЛАУ, то преимущественным оказывается использование именно этого метода. За счет меньшего времени, затрачиваемого на итерацию, более эффективным оказывается расчет нескольких дополнительных, но более экономичных с вычислительной точки

зрения итераций. Так, если для расчетной сетки рабочего колеса ступени насоса итоговые времена для SAMG(2) и SAMG(3) оказываются практически идентичными и на 25% меньшими по сравнению с САМС, то для сетки направляющего аппарата преимущество SAMG(2) более существенно по сравнению как с САМС, так и с SAMG(3).

Проведено исследование масштабируемости трех комбинаций итерационного метода подпространства Крылова с многосеточными пре-добуславливателями. Полученные результаты зависимости ускорения времени расчета одной итерации метода BiCGStab с многосеточным предобуславливателем для матрицы размером 2.7 млн. строк приведены на рис. 3. В целом графики масштабируемости дублируют друг друга, выделяя незначительное превосходство SAMG-методов по сравнению с САМС. На 40 узлах вычислительной системы «Зилант» получено ускорение в 34-37 раз, что свидетельствует о хорошей масштабируемости методов (в расчетах было использовано по одному вычислительному процессу на узел). Для второй тестовой задачи аналогичный график полностью повторяет поведение приведенных зависимостей.

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

Рис. 3. Ускорение времени расчета одной итерации метода BiCGStab с многосеточным предобуславливателем, матрица 2.7 млн. строк

Особенности реализации многосеточных методов на распределенных системах с графическими ускорителями

Выбор формата хранения данных. Для рассматриваемого набора численных методов наиболее вычислительно трудоемкой операцией оказываются матричные операции, и в частности, произведение разреженной матрицы на вектор (SpMV). В [3] был проведен предварительный анализ, который показал, что скорость выполнения данной операции на графических ускорителях значительным образом зависит от выбора формата хранения матрицы и выбор оптимального формата неоднозначен. Формат CSR, относящийся к наиболее универсальным и вместе с тем простым с точки зрения реализации на центральных процессорах, не всегда обеспечивает оптимальную производительность для этой операции на GPU. Более перспективной оказывается реализация матричных операций для формата ELLPACK, который позволяет в полной мере использовать ресурсы ускорителя. Проблема, однако, заключается в том, что ELLPACK как формат представления данных сам по себе эффективен только для узкого класса матриц с фиксированным количеством ненулевых элементов в строках. Существующие модификации этого формата позволяют несколько расширить область его применимости, но остаётся актуальным вопрос о том, как априорно определить, для каких матриц ELLPACK-формат или его модификации обеспечат лучшие результаты GPU-реализации алгоритма умножения матрицы на вектор. Замечено только, что для матриц верхних уровней многосеточного метода преимущества ELLPACK-форматов проявляются чаще и отчётливее. Таким образом, напрашиваются два вывода о форматах хранения разреженных матриц:

• использование отличных от CSR форматов хранения матриц для алгоритмов на GPU оправ-

данно, и критерии, при которых их применение имеет смысл, стоят усилий по их разработке;

• гибридные форматы хранения матриц, комбинирующие такие форматы, как ELLPACK и DIA с CSR или COO, достаточно перспективны, в первую очередь, именно для GPU-реализаций алгоритмов.

Накладные расходы в гетерогенных системах с графическими ускорителями. Запуск вычислительных ядер на графическом ускорителе и копирование данных из/в память ускорителя зачастую являются двумя основными источниками накладных расходов при использовании в расчетах графических ускорителей. Эти две операции требуют значительного времени, поскольку включают в себя целый ряд действий, предполагающих накладные расходы, такие как: взаимодействие с драйвером устройства, планирование выполнения вычислительного ядра и выполнение обмена по каналам шины PCI-Express.

Приведём пример, иллюстрирующий то, в какой степени описываемые накладные расходы могут быть существенными для рассматриваемого в работе класса задач. Обратим внимание, каким образом выполняется операция умножения разреженной матрицы на вектор в случае распределенных вычислений. В этом случае матрица хранится по процессам распределённой вычислительной системы, будучи разбитой по строкам. В соответствии с этим разбиением хранится и вектор-множитель. Для выполнения умножения алгоритму потребуются компоненты вектора-множителя, которые не являются локально хранимыми данным процессом. Процесс должен получить эти компоненты от других вычислительных MPI-процессов. Значит, именно по мере получения этих компонент и должна выполняться операция соответствующего частичного произведения матрицы на вектор. С учетом целесообразности получения нелокальных значений вектора-множителя по-

-L2: 428К “Нижеголь1

- - fl- - • L2: 428К, "Ломоносов1

» L3: 181К, Нижеголь'1

. . +- ., L3: 181К, ' Ломоносов' A L4: 75К, “Нижеголь"

- - А- - • L4: 7БК, “Ломоносов"

О 4 8 12

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

Рис. 4. График масштабируемости операции умножения разреженной матрицы на вектор на вычислительных системах с GPU-ускорителями «Нижеголь» и «Ломоносов» для матриц нескольких уровней. Рабочее колесо, 2.7 млн. неизвестных

средством МР1 большими блоками наиболее экономичным способом организации алгоритма представляется разбиение локального фрагмента матрицы данного процесса дополнительно на блоки по столбцам так, чтобы каждый блок соответствовал компоненте вектора-множителя соседних процессов (назовём такие блоки сегментами). В общем случае число таких сегментов будет равно числу задействованных в расчете МР1-процессов, хотя некоторые такие сегменты могут оказаться пустыми и не будут содержать ни одного ненулевого элемента матрицы. При увеличении количества вычислительных процессов вдвое пропорционально уменьшается объем вычислений, приходящихся на каждый процесс. Вдвое увеличивается общее количество сегментов, и на практике несколько увеличивается количество непустых сегментов локального блока матрицы. Вычисления для каждого такого сегмента приходится выполнять в отдельно запускаемом вычислительном СиОА-ядре, неизбежно увеличивая накладные расходы на запуск этих ядер, т.е. в этом случае наблюдаются одновременно два эффекта: уменьшение времени выполнения одного вычислительного ядра и увеличение количества запускаемых вычислительных ядер. Это приводит к увеличению доли накладных расходов на запуск ядер по сравнению с общим временем расчета. Исходя из такого характера накопления доли накладных расходов можно предсказать, что масштабируемость описываемой вычислительной схемы при увеличении числа процессов должна деградировать независимо от других факторов, ограничивающих масштабируемость.

Согласование копирования данных в память устройства с MPI-коммуникациями. Одной из трудностей программирования графических ускорителей можно назвать тот факт, что основные реализации библиотеки МР1 не взаимодейст-

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

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

4. Тестовые расчёты выполнены для нескольких матриц, возникающих на различных уровнях иерархии многосеточного метода CAMG для тестовой задачи размером 2.7 млн. неизвестных. Для матрицы четвёртого уровня (табл. 2) с числом неизвестных порядка 75 тысяч особенно заметен эффект деградации масштабируемости. Время выполнения операции умножения матрицы на вектор на 12 узлах с графическими ускорителями (в расчете использовалось по одной карте с узла) составляет порядка 400 микросекунд, и в процессе вычисления запускается до 10 вычислительных CUDA-ядер. Запуск ядер связан с накладными расходами, нижнюю границу для которых можно оценить следующим образом. Экспериментально, путем запуска большого количества «пустых» ядер, было получено характерное время запуска для недиагональных блоков матриц, которое составило порядка 2-3 микросекунд, и аналогичное значение для диагонального блока матрицы около 4-5 микросекунд. Суммарно это позволяет получить минимальную оценку потерь времени на запуск ядер на уровне 25-30 микросекунд. С верхней оценкой накладных расходов на запуск ядер возникают большие затруднения, поскольку такая оценка зависит от конфигурации реально запускаемого ядра, по крайней мере, в

отношении числа используемых регистров и разделяемой памяти. Накладные расходы на копирование данных оценить проще - они составляют порядка 5 микросекунд на каждую операцию копирования. Таким образом, упомянутые накладные расходы, по нашим оценкам, могут составлять, как минимум, 15% от общего времени выполнения операции для этой матрицы, и это является только оценкой «снизу».

Заключение

В статье рассмотрены вопросы применения различных многосеточных методов для ускорения сходимости итерационных методов подпространства Крылова при решении СЛАУ для эллиптических уравнений. Показано, что использование SAMG-методов в качестве предобу-славливателей, несмотря на чуть худшие показатели скорости сходимости по сравнению с CAMG, может оказаться предпочтительным за счет существенно меньшего времени расчета каждой итерации. Обсуждается ряд технологических сложностей, возникающих при реализации многосеточных методов на распределенных системах с графическими ускорителями. Показана масштабируемость параллельной реализации одного из базовых алгоритмов метода (умножения разреженной матрицы на вектор), приведены оценки накладных расходов, связанных с обменами данных и запуском вычислительных ядер, ограничивающих этот показатель.

Перспективой дальнейшей работы является исследование альтернативных подходов для агрегирования элементов во время построения матриц перехода между уровнями, оценка качества получаемых множеств элементов с точки зрения скорости сходимости итерационного процесса решения СЛАУ и разработка распределенной гибридной реализации многосеточного метода.

Требуется также дополнительная работа по поиску путей преодоления технологических

сложностей реализации базовых операций метода на GPU-ускорителях. Перспективным направлением поиска можно назвать исследование альтернативных гибридных форматов хранения матриц, используемых в ходе вычислений.

Работа частично поддержана грантом РФФИ 12-01-31002-мол-а. Расчеты выполнены на вычислительных системах «Ломоносов» Суперкомпьютерного комплекса МГУ имени М.В. Ломоносова, «Зилант» ЗАО «Т-Сервисы» и «Нижеголь» Белгородского государственного университета.

Список литературы

1. Saad Y. Iterative methods for sparse linear systems, 2-nd edition. Philadelphia: SIAM, 2003. 528 p.

2. Trottenberg U., Oosterlee C.W., Schuller A. Multigrid. N.Y.: Academic Press, 2001. 631 p.

3. Краснопольский Б.И., Медведев А.В. О решении систем линейных алгебраических уравнений на многоядерных вычислительных системах с графическими ускорителями // Параллельные вычислительные технологии (ПаВТ'2013): труды международной научной конференции (1-5 апреля 2013 г., г. Челябинск). Челябинск: Издательский центр ЮУрГУ, 2013. C. 409-420.

4. Falgout R.D. An Introduction to Algebraic Multigrid // Computing in Science and Engineering. 2006. V. 8. Is. 6. P. 24-33.

5. Vanek P., Mandel J., Brezina M. Algebraic multigrid based on smoothed aggregation for second and fourth order problems // Computing. 1996. V. 56. Is. 3. P. 179-196.

6. Van der Vorst H.A. BI-CGSTAB: A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. V. 13. № 2. P. 631-644.

7. Hypre: a library of high performance preconditioners. URL: https://computation.llnl.gov/casc/linear_ solvers/sls_hypre.html (дата обращения: 30.05.2013).

8. Bell N., Garland M. Cusp: Generic Parallel Algorithms for Sparse Matrix and Graph Computations. URL: http://cusp-library.googlecode.com (дата обращения: 30.05.2013).

9. Bell N., Dalton S., Olson L. Exposing FineGrained Parallelism in Algebraic Multigrid Methods. NVIDIA Technical Report NVR-2011 -002, 2011.

ALGORITHMIC FEATURES OF A MULTIGRID SLAE SOLVER DEVELOPMENT FOR COMPUTER SYSTEMS WITH GPU ACCELERATORS

B.I. Krasnopolsky, A V. Medvedev

A comparison of SAMG and CAMG multigrid methods to solve systems of linear algebraic equations obtained at discretization of elliptic equations is carried out on practical examples. The key peculiarities and implementation details of Krylov subspace iterative methods and multigrid methods for computer systems with GPU accelerators are discussed.

Keywords: elliptic equations, systems of linear algebraic equations, multigrid methods, graphics accelerators.

References 2. Trottenberg U., Oosterlee C.W., Schuller A. Mul-

tigrid. N.Y.: Academic Press, 2001. 631 p.

1. Saad Y. Iterative methods for sparse linear sys- 3. Krasnopol'skij B.I., Medvedev A.V. O reshenii

tems, 2-nd edition. Philadelphia: SIAM, 2003. 528 p. sistem linejnyh algebraicheskih uravnenij na mnogoya-

dernyh vychislitel'nyh sistemah s graficheskimi uskorite-lyami // Parallel'nye vychislitel'nye tekhnologii (PaVT'2013): trudy mezhdunarodnoj nauchnoj konferencii (1-5 aprelya 2013 g., g. Chelyabinsk). Chelyabinsk: Izda-tel'skij centr YuUrGU, 2013. C. 409-420.

4. Falgout R.D. An Introduction to Algebraic Multigrid // Computing in Science and Engineering. 2006. V.

8. Is. 6. P. 24-33.

5. Vanek P., Mandel J., Brezina M. Algebraic multigrid based on smoothed aggregation for second and fourth order problems // Computing. 1996. V. 56. Is. 3. P. 179-196.

6. Van der Vorst H.A. BI-CGSTAB: A fast and

smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1992. V. 13. № 2. P. 631-644.

7. Hypre: a library of high performance preconditioners. URL: https://computation.llnl.gov/casc/li-near_ solvers/sls_hypre.html (data obrashcheniya: 30.05.2013).

8. Bell N., Garland M. Cusp: Generic Parallel Algorithms for Sparse Matrix and Graph Computations. URL: http://cusp-library.googlecode.com (data obrashcheniya: 30.05.2013).

9. Bell N., Dalton S., Olson L. Exposing FineGrained Parallelism in Algebraic Multigrid Methods. NVIDIA Technical Report NVR-2011 -002, 2011.