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

Вычислительные технологии

Том 3, № 1, 1998

НЕКОТОРЫЕ ОЦЕНКИ ЭФФЕКТИВНОСТИ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ РЕШЕНИЯ СЛАУ НА ПОДПРОСТРАНСТВАХ КРЫЛОВА *

М.Ю. БАЛАНДИН, Э.П. ШУРИНА Новосибирский государственный технический университет, Россия

e-mail: shurina@fpm.nstu.nsk.su

The ABR1ORT method for solving densed ill-conditioned systems of linear equations (SLEs) and the GMRES method for solving sparsed non-symmetric SLEs are discussed; both them belong to Krylov Sequence Methods. The analysis of numerical complexity and memory requirements is performed for these methods; parallel algorythms for multiprocessor computers are briefly discussed. For these algorythms, asymptotic estimations of parallel speedup are shown.

1. Введение

Решение ряда задач математической физики (задачи гидрогазодинамики, расчета электромагнитных полей, уравнения Максвелла, Навье — Стокса и др.) методами конечных элементов (FEM) и конечных объемов (FVM) — особенно на неструктурированных сетках — приводит к системам линейных алгебраических уравнений с разреженными матрицами большой размерности, которые, как правило, являются несимметричными [10]. Частным случаем несимметричности является симметричный портрет с асимметрией коэффициентов, т.е. одинаковое расположение различных коэффициентов матрицы по ее верхнему и нижнему треугольникам [5]. В задачах гидрогазодинамики исключение давления нередко приводит к необходимости обращения близких к вырожденным (и, как следствие, плохо обусловленных) СЛАУ с плотной матрицей; плохая обусловленность может также вызываться погрешностями задания коэффициентов матрицы, связанными с неточностями измерений и/или аппроксимацией. К СЛАУ с плотной плохо обусловленной матрицей приводят также задачи, формализуемые в виде интегральных уравнений 1-го и 2-го рода [2].

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

[5, 10, 6].

В настоящей работе рассматриваются метод обобщенных минимальных невязок GMRES, эффективный для решения разреженных СЛАУ с несимметричной матрицей

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-01557.

©М.Ю. Баландин, Э.П. Шурина, 1998.

[4, 6, 8], и метод АВШОКТ из класса методов, предложенных А. А. Абрамовым для решения плотных плохо обусловленных СЛАУ [1, 3, 7]. Оба метода основываются на проектировании искомого решения на образуемое системой подпространство Крылова.

Целью работы является создание параллельных алгоритмов методов и теоретическая оценка их эффективности с точки зрения увеличения скорости решения по сравнению с последовательной реализацией; алгебраические вопросы улучшения сходимости в рамках данного исследования не рассматриваются. Для АВШОИТ и GMR.ES приводятся оценки требуемого объема памяти, вычислительной сложности и асимптотического ускорения счета для параллельных реализаций по SPMD-схеме [5].

2. Методы

Будем предполагать, что решаемая СЛАУ имеет вид Ах = Ь, где А — матрица размерности п х п, а х и Ь — векторы размерности п х 1. Пусть аг — ¿-я строка матрицы А (вектор 1 х п), а Ьг — ¿-й компонент вектора Ь. Предполагается также, что процедура предобусловливания проводится отдельно, и если она необходима, то метод применяется к уже измененной СЛАУ (так называемое "внешнее" предобусловливание), после чего найденное решение пересчитывается с учетом ранее выполненного предобусловливания. (Может использоваться, например, процедура вида Ю-1/2АЮ-1/2у = Ю-1/2Ь, где Б — диагональная матрица с главной диагональю, взятой из матрицы А.) Говоря о параллельных реализациях методов, будем обозначать доступное число процессоров как р.

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

В общем виде данный класс может быть записан как уточнение решения по формулам:

Хо = 0, Хг+1 = Хг + фг+1 ¿г+1, фг

Фг

(¿г ¿г)

г+1

а) - и,- (¿г)т, Ь!+1 = Ь) - и,Фг

и,-

(а) ¿г)

(¿г ¿г) '

¿г = (Аг )т дг, Фг = (Ьг )тдг

(1)

где вектор д является параметром, выбор которого и определяет конкретный метод из класса.

Подробное описание предложенного А. А. Абрамовым класса методов и его геометрическая интерпретация приведены в работе [7].

Метод ABR1ORT получается при дг = ег = (0, 0, ... , 1г, 0, ... , 0)т. Формулы (1) при этом могут быть упрощены, так как ¿г = (Аг)тег = а\ и фг = (Ьг)тег = Щ. Кроме того, в таком варианте метод становится чувствительным к критерию останова [3], и для повышения точности целесообразно применение специального условия, определяющего допустимость

включения очередного вектора в базис подпространства Крылова:

I (Т^к пк _ Т^к пк Пк

|(03(к+1)аз'(к) 03(к)аз'(к+1) "?(к))

Кш-2

< £,

где ] (в) —номер базовой строки на 5-й итерации.

Метод GMR.ES [4, 6, 9] был предложен для решения несимметричных СЛАУ и основан на построении базиса в соответствующем системе подпространстве Крылова К3(А, Го). Затем решение уточняется некоторой добавкой, представленной в виде разложения по этому базису.

В несимметричном случае основная сложность заключается в том, что для построения ортогонального базиса не могут быть использованы рекуррентные соотношения с малой глубиной зависимости типа тех, что применяются в методе сопряженных градиентов. Поэтому в GMRES ортонормальный базис (г?} строится на основе следующего соотношения, выражающего полную зависимость всех векторов:

' 1 ^+1,^+1 = Лг» - ^ ^г?, г1 = -—- Го. (2)

3=1 -Го-

В матрично-векторных обозначениях соотношение (2) может быть записано как

ЛУг = УгИг + = УтИг, (3)

где Vi = [г1|г2|...|г^], а Н —матрица размерности (г + 1) х г в верхней форме Хессенберга.

Особенностью GMRES является то, что с целью экономии памяти выбирается некоторая размерность подпространства Крылова т (причем т ^ п), и после того как т базисных векторов построены решение уточняется. Построение очередного вектора vi обычно называют GMRES-итерацией, а т таких итераций с уточнением решения — GMRES-циклом. Один цикл требует т +1 матрично-векторное умножение.

Система (3) является плотной и имеет малые размеры; кроме того, она переопределена, а потому должна решаться в смысле наименьших квадратов.

Следует также отметить, что при т = п GMRES фактически становится методом абрамовского класса. Подробное описание GMRES приведено в [6, 8, 9].

3. Анализ методов, параллельные алгоритмы и асимптотические оценки

Схема метода ABR1ORT характерна тем, что объем вычислительной работы от итерации к итерации уменьшается, так как в качестве векторов ^ последовательно выступают строки матрицы Л, которые в дальнейших расчетах уже не участвуют. Основная тяжесть вычислительной работы приходится на вычисления скалярных произведений строк матрицы и их линейное комбинирование. В общей сложности на г-й итерации вычисляется п — г + 1 скалярное произведение (г-я строка на себя и на каждую из п — г нижележащих) и производится п - г операций линейного комбинирования строк.

Несложные подсчеты показывают, что на г-й итерации с учетом проверки значимости строк выполняется в общей сложности 4п2 — 4(п — 1)г + 7п +1 арифметическое действие.

Поскольку

п

^(4п2 - 4(п - 1)г + 7п + 1) = 2п3 + 7п2 + 3п = ^А(п),

г=1

то можно говорить, что метод имеет сложность 0(п3). Действительно, его вычислительная схема обработки строк матрицы схожа со схемой метода Гаусса, который имеет такую же сложность.

Интерес представляет также исследование метода с точки зрения требуемой им машинной памяти. При программной реализации метод требует хранения матрицы СЛАУ (использование разреженных форматов нецелесообразно, так как матрица постоянно пе-ресчитывается), т.е. п2 вещественных чисел; векторов решения и правой части, т.е. 2п вещественных чисел; значений скалярных произведений фиксированной строки матрицы на саму себя и на все остальные строки, т. е. п вещественных чисел; вектора флагов значимости строк матрицы, т. е. п логических величин. Таким образом, всего требуется хранение п2 + 3п вещественных и п логических величин (объем памяти под счетчики циклов и временные переменные скалярного типа можно считать пренебрежимо малым). Если вещественное число требует 4 байт, а логическая величина 1 байт, то для решения СЛАУ размерности п х п метод ABR1ORT требует 4п2 + 13п байт; при использовании вещественных чисел двойной точности требования возрастают до 8п2 + 25п байт. Следовательно, требования метода ABR1ORT к объему доступной памяти составляют 0(п2).

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

Если решение малой переопределенной СЛАУ производится с помощью QR-факториза-ции с применением вращений Гивенса [8-10], то подсчет показывает, что общее количество арифметических действий, выполняемых на одном GMRES-цикле,

N с(п, т) = 5п + 9т + 6к + 8пт + 2пт2 + 2т2 + 2тк + 1.

С учетом того, что т ^ п (возможности современных вычислительных систем обычно ограничивают выбор размерности подпространства Крылова величиной т = 0.01п) и к = О(п) (обычно к < 100п), можно говорить, что один GMRES-цикл имеет сложность 0(пт2). Проблема оптимального автоматического выбора размерности подпространства Крылова т в настоящее время не решена, однако число GMRES-циклов теоретически не превосходит [п/т].

При выборе т = п для решения СЛАУ становится достаточно только одного GMRES-цикла, и метод приобретает сложность 0(п • п2) = 0(п3). Действительно, данная величина совпадает со сложностью метода ABR1ORT.

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

Вычислительная схема требует работы со следующими объектами: разреженной матрицей исходной СЛАУ — п + 2к вещественных и п + к +1 целых (причем обычно двойной

длины!) чисел; векторами правой части, искомого решения и невязки — 3п вещественных чисел; векторами базиса подпространства Крылова — тп вещественных чисел; плотной матрицей малой переопределенной СЛАУ — т(т +1) вещественных чисел; коэффициентами вращений Гивенса — 2т вещественных чисел; вектором решения малой СЛАУ (как с коэффициентами линейного комбинирования векторов базиса подпространства Крылова) — т вещественных чисел.

Полагая, как и прежде, объем памяти для счетчиков циклов и временных переменных скалярного типа пренебрежимо малым, можно получить, что всего требуется хранить п + к + 1 целых чисел двойной длины и 4п + 2к + тп + (т + 1)2 + 3т +1 — вещественных. При использовании обычной точности это требует 20п+12к+4тп+4(т+1)2 + 12т+8 байт, двойной — 36п + 20к + 8тп + 8(т + 1)2 + 24т + 16 байт. С учетом того, что 1 < т ^ п и к ^ п, можно говорить, что требования метода GMRES к объему доступной памяти составляют О(тп).

Как было показано ранее, основные вычислительные затраты в ABR1ORT определяются операциями нахождения скалярных произведений и линейного комбинирования век-

31

торов. Можно подсчитать, что эти действия требуют 2п3 + ^п2 — ^п арифметических операций. Доля соответствующих вычислительных затрат в схеме метода (требующей NА(п) операций), таким образом, равна

п 3 3 2 1

2п3 +— п--п

2 2

2п3 + 7п2 + 3п

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

\ п(3п2 + 2п + пр + р + 1)

^ (п) = ^--.

р

Таким образом, асимптотическое ускорение для ABR1ORT равно

т NА (п) 2

11т = -р. (4)

^(п) 3

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

Численные эксперименты (см. раздел 4) на матрицах Гильберта (наихудший случай плохо обусловленных матриц) показали результаты, довольно хорошо согласующиеся с

(4).

Вычислительная схема GMR.ES значительно более сложна. Как можно видеть из сказанного ранее, в ней присутствуют два интенсивных в вычислительном плане этапа: умножение матрицы СЛАУ на вектор (для нахождения невязок) и нахождение скалярных произведений векторов (при построении ортонормального базиса подпространства Крылова). Кроме того, именно на хранение векторов базиса {у^} в GMRES приходится большая часть затрат памяти (ши вещественных чисел, что часто оказывается больше объема памяти под разреженную матрицу СЛАУ). При этом в GMRES-схеме присутствует достаточно самостоятельная подзадача — решение переопределенной СЛАУ малой размерности в смысле наименьших квадратов, которая должна рассматриваться отдельно.

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

Используемая топология в значительной степени определяется тем, каким способом решается малая переопределенная СЛАУ. Если данная подзадача решается последовательно (что вполне возможно из-за ее малых размеров), то наиболее подходящей оказывается топология "звезда".

Пусть малая СЛАУ решается последовательно с применением метода вращений Ги-венса. Тогда анализ для топологии "звезда" показывает, что за один цикл вычислительная нагрузка центрального (наиболее занятого) процессора "звезды" составляет

с 4и + 8к + 1 + 6ши + 2шк + 2ш2и + 2ш + 2ш2р + 2шр + 2шир

"о

' р

арифметических операций.

Таким образом, асимптотическое ускорение, равное

Nс(и,ш) 1 (5 + 8ш + 2ш2)р . .

11т —РТ7-т = - • -¡5-, (5)

н^те "р (и, ш) 2 2 + 3ш + ш2 + шр

оказывается функцией не только числа процессоров, но и размерности подпространства Крылова ш (ш ^ и). Так, можно видеть, что при больших ш (ш ^ и) асимптотическое ускорение составит

N с(и,ш)'

11т 11т —^-г = Р-

н^те ут^н "р (и,ш) )

Формула (5) показывает, что выбор глубины пространства Крылова является весьма важной проблемой. Теоретически выбор большего значения должен приводить к более быстрому достижению решения, однако это справедливо лишь для вычислений в точной арифметике. Как показали численные эксперименты, при чрезмерном увеличении ш возможно даже замедление расчетов, так как метод становится чрезвычайно чувствительным

к вспомогательному критерию точности, который на заключительных циклах перестает распознавать "избыточные" векторы в базисе подпространства Крылова.

4. Некоторые результаты численных экспериментов

В табл. 1 и 2 показаны значения ускорения при параллельном счете методами ABR1ORT и GMRES в сравнении с теоретическими оценками. Эксперименты проводились на вычислительной системе PARSYTEC в среде Power Xplorer (см. также [3]). Обозначения: tp — время решения на процессорах; s™0^ — значение теоретической оценки ускорения (см. раздел 3); ЗПрак = t\/tp — практически достигнутое ускорение. Погрешность измерения времени при расчетах составляла ±1с.

Метод ABR1ORT исследовался на симметричных матрицах Гильберта, GMRES — на матрице размерности 8112, полученной из задачи расчета электромагнитного поля методом конечных элементов (подробности см. в [8]). При этом расчеты производились без пре-добусловливания и с диагональным предобусловливанием вида D-1/2AD-1/2y = D-1/2b.

Таблица 1. Метод ABR1ORT. Таблица 2. Метод GMRES.

n Р tp „теор sp „прак sp

700 1 6" 1.0 1.0

2 6" 1.3 1.0

3 5" 2.0 1.2

4 3" 2.7 2.0

800 1 7" 1.0 1.0

2 6" 1.3 1.2

3 5" 2.0 1.4

4 4" 2.7 1.8

900 1 10" 1.0 1.0

2 8" 1.3 1.3

3 6" 2.0 1.7

4 5" 2.7 2.0

m Р tp „теор „p „прак sp

б/предоб. предоб. б/предоб. предоб.

10 1 11/36// 10// 1.00 1.00 1.0

2 8'27" 7// 1.88 1.37 1.4

4 5/6// 4// 3.31 2.28 2.5

20 1 7/8// 9// 1.00 1.00 1.0

2 4/59// 6// 1.92 1.43 1.3

4 3/22// 4// 3.56 2.12 2.3

30 1 5/28// 8// 1.00 1.00 1.0

2 4/1// 6// 1.94 1.36 1.3

4 2/37// 3// 3.68 2.09 2.7

5. Заключение

В процессе исследований был детально проанализирован и эффективно реализован метод ABR1ORT, пригодный для решения плотных плохо обусловленных СЛАУ с высокой точностью (в том числе недо- и переопределенных систем), а также реализован и исследован метод GMRES; в связи со сложностью схемы GMRES представляются целесообразными эксперименты с реализацией метода на системах с разделяемой памятью [5].

Численные эксперименты показали достаточно хорошую (с учетом сравнительно низкой скорости передачи данных по межпроцессорным соединениям у вычислительной системы PARSYTEC [6]) согласованность с теоретическими оценками (4), (5)

Дальнейшие исследования в данном направлении связаны в первую очередь с применением ABR1ORT в схеме GMRES и проблемой автоматического выбора размерности подпространства Крылова, а также с использованием этих методов для решения нелинейных систем.

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

[1] Абрамов А. А. Об одном методе решения плохо обусловленных систем линейных алгебраических уравнений. Журн. вычислит. матем. и матем. физики, 31, №4, 1991, 483-492.

[2] Абрамов А. А. О свойствах метода Крейга при решении линейных некорректных задач. Там же, 35, №1, 1995, 144-150.

[3] Баландин М.Ю., Шурина Э.П., Чернышев О. В. Анализ параллельных алгоритмов решения плохо обусловленных СЛАУ. Тр. III междунар. науч.-техн. конф. "Актуальные проблемы электронного приборостроения" (АПЭП-96), Т. 6, Ч. 1, Новосибирск, 1996, 74-77.

[4] Ильин В. П. Методы неполной факторизации для решения алгебраических систем. Физматлит, М., 1995.

[5] Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. Пер. с англ. Мир, М., 1991.

[6] Aspects of Computational Science. Stichting Nationale Computer Faciliteiten, The Netherlands, 1995.

[7] BALANDIN M. Yu. ET AL. On parallelization of new algorithm for solving systems of linear equations. In "Joint Bulletin of the Novosibirsk Computing Center and the Institute of Informatics Systems", Ser. "Computer Science", 1996, 17-26.

[8] Balandin M., Chernyshev O., Shurina E. Analysis of methods for solving large-scale non-symmetric linear systems with sparsed matrices. In "Parallel Computing Technologies, 4th Int. Conf. (PaCT-97) Proc." Springer-Verlag, "Lecture Notes in Computer Science" series, 1277, 1997, 336-343.

[9] KADioGLu M., Mudrick S. On the implementation of the GMRES(m) method to elliptic equations in meteorology. J. of Comput. Phys., 102, 1992, 348-359.

[10] Natarajan R. Finite element applications on a shared-memory multiprocessor: algorithms and experimental results. Ibid., 94, 1991, 352-381.

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