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

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ...

УДК 517.968:519.612:004.021

Е. Н. Акимова, Д. В. Белоусов

ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ РЕШЕНИЯ СЛАУ С БЛОЧНО-ТРЕХДИАГОНАЛЬНЫМИ МАТРИЦАМИ НА МНОГОПРОЦЕССОРНЫХ ВЫЧИСЛИТЕЛЯХ

Для решения СЛАУ с блочно-трехдиагональными матрицами предложены и численно реализованы на многопроцессорных вычислителях различного типа параллельный алгоритм матричной прогонки и параллельный метод сопряженных градиентов с предобуславливателем. Проведено сравнение времени счета параллельных алгоритмов на видеоускорителе NVIDIA GeForce, 4-ядерном процессоре Intel Core I5-750 и многопроцессорном вычислительном комплексе МВС-1000 при решении задачи о нахождении распределения потенциала на поверхности земли в проводящей среде. Параллельные алгоритмы; эффективность и ускорение; суперкомпьютер; видеоускоритель

ВВЕДЕНИЕ

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

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

дС/ 1 д

дт г4 дг

( N-1

q

_ ~ ЭС. ^ r У D..—J-^ j Эг

j_1

(1)

У

где N - число компонентов сплава; г - пространственная координата; Ск - концентрация ]-го компонента; Оі}. - парциальные коэффициенты взаимной диффузии; q - показатель степени, зависящий от симметрии задачи: 0 - для плоской, 1 - для цилиндрической, 2 - для сферической.

При использовании абсолютно устойчивой неявной разностной схемы система дифференциальных уравнений (1) сводится к системе уравнений с блочно-трехдиагональной матри-

Контактная информация: (343)375-34-46 Статья рекомендована к публикации программным комитетом Международной научной конференции "Параллельные вычислительные технологии 2011". Работа выполнена в рамках Программы фундаментальных исследований Президиума РАН № 14 при поддержке УрО РАН, проект 09-П-1-1003.

цей [1-2]. Для /-го компонента в к-м узле пространственной сетки система уравнений имеет вид

N-1

Х[ (к-1)+(к) +Ь£: (к+1)]=<, (2)

]=1

где а/к, Ьк, с/к, dгk - коэффициенты, рассчитываемые в зависимости от модели, Ск(к) - концентрация к-го компонента в к-м узле сетки.

Важнейшими задачами исследования неоднородности земной коры являются задачи электроразведки. Одним из известных методов электроразведки является метод вертикального электрического зондирования (ВЭЗ). На поверхности земли собирают электроразведочную установку, которая состоит из двух питающих и двух приемных электродов. К питающим электродам А и В подключают источник тока, а на приемных электродах М и N измеряют напряженность электрического поля. По результатам выполненных измерений вычисляют кажущееся электрическое сопротивление рк = = КЬУЖ / 1АВ для неоднородной среды, где К -геометрический коэффициент, зависящий от расстояний между электродами А, В, Ми N, ЬУ^ - разность потенциалов на приемных

электродах М и N, 1АВ - сила тока, протекающего в питающей линии.

Для модели среды с плоско-параллельными горизонтальными границами, когда слои земной коры залегают горизонтально и сопротивление среды зависит только от глубины р = р(г), кажущееся сопротивление р = р(г), где г = АО -расстояние от центра приемной цепи до питающего электрода. Задача интерпретации результатов ВЭЗ заключается в определении р(г), дающего «электрический разрез» среды по из-

вестным значениям р(г) (ось г направлена вниз, ось г направлена вправо).

Задача ВЭЗ о нахождении потенциала У0 (г, г = 0) на поверхности земли сводится к решению двумерного уравнения Лапласа в цилиндрической системе координат [3]

АТ. Э2К 1ЭК Э 2V п

АК °—— :------------:—-_ 0

Эг 2 r Эг Эz

(3)

с граничными условиями непрерывности потенциала и непрерывности нормальной составляющей к границам плотности тока:

V

z _ І _V1

1

эк

z _ Г'

r0 Эz

1

эк

z _ І

z _ І

условия на горизонтальных прямых z = І, ЭК , ЭК ,

—0 z_0 _ 0 - сверху, —0 r_0 _ 0 - слева, К0 = 0 -

Эz

справа и снизу. Здесь

Эг

p _-

2nr2

ЭК

Э;

(4)

После использования конечно-разностной аппроксимации краевая задача (3)-(4) сводится к решению СЛАУ с блочно-трехдиагональной матрицей (рис. 1).

Другой важной задачей электроразведки является задача бокового каротажного зондирования (БКЗ), предусматривающая использование приборов однотипных зондов разной длины при измерении потенциала электрического поля. В результате интерпретации данных каротажа получают значение удельного электрического сопротивления пласта, близкое к истинному.

В работе [4] показано, что после использования конечно-разностной аппроксимации задача БКЗ сводится к решению системы линейных алгебраических уравнений с блочнотрехдиагональной матрицей следующего вида:

raV I -

(bVz )z

_ F,

(5)

где А - матрица размерности N х 2Nz) х N х х 2Nz), имеющая блоки размерности Nг х Nг (рис. 1); У и ^ - векторы размерности Nг х 2Nz; Nг - число точек по г; 2 х Nz - число точек по г; о = 1 / р - коэффициент электропроводности;

а.. _ о

j

h

(r)

h

(z ) Л

; --

b.. _ о

ij

; :-

,z.:-

j

(z)

, z. --

j

(К );,„ _(к - v_u)/ a; ;

(v);,„ _ (v>l,j - к,. ) / й ‘;>.

(К),„ - Vj-l)/a'z

(v = (К/,.» - V,.)/ й'Д

й(') _(a<;):A«)/2, Л' _; -r-l, i_ 1,..

йf _(j :j)/2, j _z, -z,_l,

, N.

j _-Nz:1,...,Nz.

Рис. 1. Вид матрицы СЛАУ

Целью данной работы является разработка и реализация алгоритмов решения СЛАУ с блочно-трехдиагональными матрицами на многопроцессорных вычислителях различного типа: на видеоускорителе NVIDIA GeForce

GTX 285, 4-ядерном процессоре Intel Core I5-750 и многопроцессорном вычислительном комплексе МВС-1000 и сравнение времени счета параллельных алгоритмов при решении модельной задачи.

1. ПАРАЛЛЕЛЬНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ

Рассмотрим систему уравнений с блочнотрехдиагональными матрицами

'од-ад = F,, i=0,

«-4Y-1 + CY-в?1+1 = Ft, i=1,...,N-1, (6)

-ANYN-1 + CnYn = Fn , i = N,

где Yi - искомые векторы размерности n, Fi -

заданные векторы размерности n, Ai, Bi, Ci -квадратные матрицы порядка n.

Для решения СЛАУ (6) предлагается использовать параллельный прямой метод матричной прогонки [5] и параллельный итерационный метод сопряженных градиентов с предо-буславливателем в случае решения СЛАУ

с симметричной положительно-определенной

матрицей.

o

r

r

2

2

Ранее параллельный метод матричной прогонки, реализованный на многопроцессорном вычислительном комплексе МВС-1000, эффективно использовался при решении модельной задачи диффузии о насыщении плоской металлической пластины тремя компонентами с различной диффузионной подвижностью, а также при моделировании эволюции выделений в двухфазном многокомпонентном сплаве [6-7].

1.1. Прямой метод решения задачи

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

При построении параллельного алгоритма исходную область Р разобьем на Ь подобластей (рис. 2) вертикальными линиями так, что N = Ь х М.

Рис. 2. Разбиение области на подобласти

В качестве параметрических неизвестных выберем векторы УК , К = 0, М,....А, связывающие неизвестные на сетке по вертикали. Относительно УК строится редуцированная система уравнений. Для этого в подобластях, определяемых интервалами (К, К + М), рассматриваются задачи:

-4Ц-1+СЩ - цим = о, ик =

и1 =

^к+м

_ _ _ _ о

-и +си - ви+1 = о, ик =о

її

(7)

ип = ■■ •

^к+м о ’

-4У-+СУ - ял+1 = о, У =

о Vі =

■■ 5 ¥ к+м

< \ о Ґ \ о

В.У+= о, уп = +1 к о 1о, Уп = к+м о I1,

(8)

-4^-1 +сд-цщ+1 = р, ш =

ш =

"к+м

(9)

где і = к + 1,..., к + м- 1.

Утверждение 1. Если ии? - решения

задач (7), V',■■■,V? - решения задач (8), Ж -

решения задачи (9), а Уі - решения исходной зада чи (6) на (к, к + м), тогда

У=(иий у +(уу ‘..у; )Гкм + ш (1о)

После подстановки соотношений (Ю) в исходную систему (6) в точках к = о, м,...,Ы, получим систему уравнений относительно параметров Ук. Эта система уравнений по своей структуре аналогична (6), имеет меньшую размерность и следующий вид:

ЦСо -Вои]Уо-[ВоУ ]Ум = Р + В,% к = о; -[4ик-1 ]Ук-м +[Ск -4кУк-1 -Вкик+1 ]Ук -

-[ ВкУк+1 ]ук + м = Рк + 4кшк-1 + ВкЖк +1,

к = м ,2м ,■■■, N - м;

-[4Л-1 У-м + С -4мум-1 У =

= Р + -1, к = N,

(11)

где ик и Ук - квадратные матрицы порядка п.

Задача (11) решается классическим алгоритмом матричной прогонки [8] на одном процессоре, задачи (7)-(9) решаются независимо в Ь подобластях на Ь процессорах.

Матрицы иі, Уі и вектор Жі на интервалах (к, к + м) вычисляются независимо на Ь процессорах по следующим формулам.

а) Прямой ход прогонки:

Р к +1 = СА'+14к+1,

а к+1 Ск+1Вк+1,

У к +1 Ск+1Рк +1, аі [Сі 4а-1] Ві

Рі = [С - 4 а-1 ]-14 Рі-1,

їі = [Сі- 4 а-1 ]-1 [р+4 у--1 ], і = к + 2,..., к + м -1.

б) Обратный ход прогонки:

ик +м -1 = Р

(12)

к +м -1>

Ук +м-1 ак +м -1,

ЩК +М -1 = Ук +М -1, и 1 = аи 1+1 + Рг, У = аг'У +1,

щ = аг^+1 + Уг, г = К + М - 2,..., К +1.

(13)

После вычисления параметров УК остальные искомые неизвестные находятся по формуле

о

о

о

о

о

о

п

(10) также независимо на каждом из Ь интервалов на Ь процессорах.

Схема параллельного алгоритма имеет вид:

(12) ^ (13) ^ (11) ^ (10).

Утверждение 2. Если для исходной системы (6) выполняются достаточные условия устойчивости метода матричной прогонки по А. А. Самарскому [8]

|с0-ч| < 1, са4 < 1,

ЦСдЦ+||сг1д|| < 1, г=1,..., N -1,

причем хотя бы одно из неравенств - строгое, то эти же условия достаточны и для устойчивости метода матричной прогонки при решении системы уравнений (11) относительно параметров УК (см. [5]).

1.2. Итерационный метод решения задачи

Одним из быстрых итерационных методов решения СЛАУ с симметричной положительноопределенной матрицей является метод сопряженных градиентов (МСГ) [9].

Введение предобуславливания применяется с целью ускорения сходимости итерационного процесса и состоит в том, что исходная СЛАУ

Ах = Ь (14)

заменяется на систему уравнений

С~хАх = СХЬ, (15)

для которой итерационный метод (в нашем случае МСГ) сходится существенно быстрее.

Пусть С = С > 0. Предположим, что предо-буславливатель представлен в виде С = В В, где В - невырожденная квадратная матрица. Умножим слева обе части системы (15) на В и положим у = Вх. В результате придем к эквивалентной системе уравнений

Ау=Ь, (16)

где А = (В*)-1 АВЬ = (В*)-1 Ь.

Условием выбора предобуславливателя С является следующее:

еопё( А) << еопё( А),

еопа(л)=Ь2^, еопа(4)=^

(17)

где еопё(А) и еопё( А) - числа обусловленности

матриц А и А; X1тах и XШ1п, хтт - наибОльшее и наименьшее собственные значения матриц А и А, соответственно.

Для системы уравнений (15) метод сопряженных градиентов с предобуславливателем С имеет следующий вид [10]:

г0 = Ь - Ах0, р° = С~хт0, / = р0,

(гк, / )

хк+1 = хк + а крк,

а

к (4рк, рк)

(18)

гк+1 = гк - ак4рк, -к+1 = С~1гк+1,

рк+1 -к+1 , в рк « (гк+1, -™)

Р =- + РкР , Рк = ( к ГЛ •

(г , - )

Условием останова итерационного процесса

является

\4х - Ь11 И

< 8.

(19)

Распараллеливание итерационного метода сопряженных градиентов основано на разбиении матрицы 4 горизонтальными полосами на т блоков, а вектора решения х и вектора правой части Ь СЛАУ на т частей так, что п = т х Ь, где п - размерность системы уравнений, т -число процессоров, Ь - число строк матрицы в блоке (рис. 3).

Рис. 3. Схема распределения данных по процессорам

На каждой итерации каждый из т процессоров вычисляет свою часть вектора решения. В случае умножения матрицы 4 на вектор х каждый из т процессоров умножает свою часть строк матрицы 4 на вектор х. Но8І;-процессор отвечает за пересылки данных и также вычисляет свою часть вектора решения.

2. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

Параллельный алгоритм матричной прогонки (ПАМП) и параллельный метод сопряженных градиентов с предобуславливателем (ПМСГ) реализованы на следующих современных многопроцессорных вычислительных системах: многопроцессорном вычислительном

|| y^ — Y "пЦ

уы\

(20)

комплексе МВС-1000/64 с помощью технологии MPI [11], 4-ядерном процессоре Intel Core I5-750 (CPU) на языке С++ в среде разработки «Visual Studw» (распараллеливание по потокам данных) и видеоускорителе GeForce GTX 2В5 (GPU) с помощью технологии CUDA [12].

Многопроцессорный вычислительный комплекс МВС-1000 - российский массивно-

параллельный суперкомпьютер кластерного типа с распределенной памятью, установленный в Институте математики и механики УрО РАН (г. Екатеринбург).

МВС-1000/64 состоит из 14 2-процессорных 2-ядерных модулей AMD Opteron 64 bit (2.6 Ггц), интерфейса GbitEthernet и 112 Гб оперативной памяти.

Технические характеристики вычислительной платформы Intel Core I5-750 с видеоускорителем NVIDIA GeForce GTX 2В5 приводятся в табл. 1.

Таблица 1

Технические характеристики

вычислительной плат ы м р о

4-ядер. Intel

СРи Core I5-750

Частота проц. (ГГ ц) 2.66

Оператив. память (Гб) В

Разрядность ОС (Бит) 64

GeForce

ОРИ GTX 285

Количество проц. ядер 240

Частота ядра (МГ ц) 64В

Частота проц. (МГц) 1476

Колич. видеопамяти (Мб) 1024

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

Исходные данные и модельное решение задачи предоставлены лабораторией скважинной геофизики Института нефтегазовой геологии и геофизики СО РАН (г. Новосибирск).

После дискретизации задача сводится к решению СЛАУ с плохо обусловленной симметричной положительно-определенной блочнотрехдиагональной матрицей вида (6) размерности 76136x76136 с квадратными блоками порядка 248 (см. рис. 1).

Приближенное решение задачи сравнивалось с модельным решением с помощью вычисления относительной погрешности

где Ум - модельное решение задачи, УП - приближенное решение задачи.

и-потенциал (решение СЛАУ)

О 100 200 300 400 500 600 700 800 900 1000

Рис. 4. Решение модельной задачи

Условие (20) выбиралось в качестве критерия останова итерационного ПМСГ при решении модельной задачи.

Предварительно находилось число обусловленности исходной матрицы А

сопа(А) = ^тах / 1,3-10",

^тах = 1,4-106, ^тт = 1,1 ■ 10-5 > 0,

где ^тах и ^тт - наибольшее и наименьшее собственные значения исходной матрицы.

В случае решения задачи предобусловлен-ным ПМСГ с целью проверки условия (17) находилось число обусловленности матрицы А

соп^А) = » 4,1 • 109 < соп^А).

^ тт

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

Бт = Т / Гт, Ет = Бт / т, Б = Т / Тг, где Тт - время выполнения параллельного алгоритма на МВС-1000 либо на многоядерном процессоре с числом процессоров или ядер т (т > 1), Т - время выполнения последовательного алгоритма на одном процессоре либо на одном ядре, Т2 - время решения задачи на видеоускорителе. Тт представляет собой совокупность чистого времени счета и накладных расходов на межпроцессорные обмены Тт = Тс + Та. Число процессоров т соответствует упомянутому разбиению векторов на т частей и разбиению исходной области на т подобластей.

В общем случае эффективность распараллеливания меняется в пределах 0 < Em < 1. В идеальном случае при равномерной и сбалансированной загрузке процессоров и минимальном времени обменов между ними Em близко к единице, но при решении практических задач она уменьшается за счет накладных расходов.

В табл. 2 и 3 приведены времена счета и коэффициенты ускорения и эффективности решения модельной задачи на 4-ядерном процессоре Intel Core I5-750, видеоускорителе NVIDIA GeForce GTX 285 и многопроцессорном комплексе МВС-1000/64 с помощью параллельного метода сопряженных градиентов с предобуслав-ливателем и параллельного алгоритма матричной прогонки при оПАМп ~ 2-10-7. Использование технологии OpenMP для 4-ядерного процессора Intel Core I5-750 с общей памятью уменьшает время решения задачи.

Т аблица 2 Решение задачи методом ПМСГ

Вычислитель Tm, с Sm (S) E J-,m

Intel Core I5-750 (1 ядро) 57 21-OpenMP — —

Intel Core I5-750 (2 ядра) 46 16-OpenMP 1,24 0,62

Intel Core I5-750 (4 ядра) 36 14-OpenMP 1,50 0,40

GeForce GTX 285 16 3,56 —

МВС—1000/64 (1 проц.) 120 — —

МВС—1000/64 (2 проц.) 65 1,85 0,92

МВС—1000/64 (4 проц.) 34 3,53 0,88

Т аблица 3 Решение задачи методом ПАМП

Вычислитель Tm, с Sm (S) E J-,m

Intel Core I5-750 (1 ядро) 52 21-OpenMP — —

Intel Core I5-750 (2 ядра) 28 18-OpenMP 1,86 0,93

Intel Core I5-750 (4 ядра) 16 3,25 0,81

GeForce GTX 285 10 5,2 —

МВС—1000/64 (1 проц.) 96 — —

МВС—1000/64 (2 проц.) 60 1,6 0,80

МВС—1000/64 (4 проц.) 31 3,1 0,77

Заметим, что время решения задачи с помощью метода сопряженных градиентов без пре-добуславливателя на одном ядре Intel Core I5-

750 при оМСГ = 10-3 составило 55 минут, что существенно превышает времена решения задачи, представленные в табл. 2 и табл. 3.

ЗАКЛЮЧЕНИЕ

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

Проведено сравнение времени счета параллельных алгоритмов на видеоускорителе NVIDIA GeForce GTX 285, 4-ядерном процессоре Intel Core I5-750 и многопроцессорном вычислительном комплексе МВС-1000/64 при решении модельной задачи о нахождении распределения потенциала на поверхности земли в проводящей среде.

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

Авторы выражают признательность за полезные советы и внимание к работе члену-корреспонденту РАН Владимиру Васильевичу Васину.

СПИСОК ЛИТЕРАТУРЫ

1. Crank J. The Mathematics of Diffusion. Oxford: Clarendon press, 1975.

2. Самарский А. А. Теория разностных схем. М.: Наука, 1983.

3. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1966.

4. Дашевский Ю. А., Суродина И. В., Эпов М. И. Квазитрехмерное математическое моделирование диаграмм неосесимметричных зондов постоянного тока в анизотропных разрезах // Ch6. ЖИМ. 2002. Т. 5. №3 (11). С. 76-91.

5. Акимова Е. Н. Распараллеливание алгоритма матричной прогонки // Математическое моделирование. 1994. Т. 6. № 9. С. 61-67.

6. Акимова Е. Н., Горбачев И. И., Попов В. В. Решение задач многокомпонентной диффузии с помощью алгоритма матричной прогонки // Мат. моделирование. 2005.Т. 17. № 9. C. 85-92.

7. Горбачев И. И., Попов В. В., Акимова Е. Н. Моделирование диффузионного взаимодействия карбонитридных выделений с аустенитной матрицей с учетом возможности изменения их состава // Физика металлов и металловедение. 2006. Т. 102. № 1. С. 22-32.

8. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

9. Фаддеев В. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. М.: Гос. издат. физ.-мат. литературы, 1963.

10. Амосов А. А. Циркулянтно предобуслов-ленный метод сопряженных градиентов и его применение для численного решения интегрального уравнения переноса излучения // Современные проблемы математического моделирования: Тр. XI Всерос-сийск. шк.-семинара. Ростов-на-Дону: Изд-во Рос-товск. гос. ун-та, 2005. Выпуск 4. С. 49-65.

11. Руководство пользователя системы МВС-1000 / А. В. Баранов [и др.]. ИКЬ: Ы1р://рага1-1e1.ru/mvs/user.htm1 (дата обращения: 15.02.2011).

12. Берилло А. КУГО1А СИБА - неграфические вычисления на графических процессорах. ИКЬ: Ы1р://'№№^1хЫ.сотМёео3/сиёа-1^Мт1 (дата обращения: 15.02.2011).

ОБ АВТОРАХ

Акимова Елена Николаевна, вед. науч. сотр. отдела некорректных задач анализа и приложений Ин-та математики и механики УрО РАН, проф. каф. вычислительн. методов и уравнений матем. физики Уральск. федерального ун-та, г. Екатеринбург. Дипл. математик (Уральск. гос. ун-т, г. Екатеринбург, 1982). Д-р физ.-мат. наук по матем. моделированию, численным методам и комплексам программ (ЮУрГУ, 2009). Иссл. в обл. параллельных алгоритмов решения прикладных задач.

Белоусов Дмитрий Владимирович, вед. программист отдела вычислительных сетей Ин-та математики и механики УрО РАН, г. Екатеринбург. Дипл. матем.-программист по матем. обеспечению и администрированию инф. систем (УГТУ-УПИ, 2006). Иссл. в обл. параллельных вычислений на многопроцессорных системах.