авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 || 3 | 4 |

«Серия Суперкомпьютерное Образование Координационный совет Системы научно-образовательных центров суперкомпьютерных технологий ...»

-- [ Страница 2 ] --

Простой блочный алгоритм Распределение данных между процессорами производится сле дующим образом. Матрица A распределяется по процессорам бло ками, содержащими (n / p) n элементов, B – n ( n / p ) элементов.

Далее проводится матричное умножение соответствующих блоков, и определяются блоки матрицы C размером ( n / p) (n / p), стоящие на главной диагонали.

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

A1 B1 B2 B3 C A2 C A3 C шаг С A3 B1 B2 B С A С A шаг A2 B1 B 2 B3 C A3 C A1 C шаг Заметим, что для выполнения этого алгоритма оперативная па мять каждого ПЭ должна вмещать не менее 3n2 / p чисел с плаваю щей точкой. Тогда общее количество пересылаемых элементов мат рицы A :

n p (этапов) (размер блока Ai )= n2.

p Предполагается, что пересылка блоков Ai, i 1,..., p между процес сорными элементами на каждом шаге осуществляется одновременно.

Алгоритм Кэннона Для этого алгоритма количество процессоров p s 2, где s – натуральное число. Матрицы A и B делятся на p блоков согласо ванных размеров: Aij, Bij, i, j 1,..., p. На каждый процессорный элемент pij рассылаются соответствующие матрицы Aij, Bij. Нахо дится их произведение. Затем на каждой итерации производится: пе ресылка матрицы Aij левому соседу, а матрицы Bij – верхнему в двумерной решетке процессоров, в которой первый и последний уз лы как в ряду решетки, так и в ее колонке имеют прямое соединение;

вычисление произведения полученных матриц. После выполнения p итераций получается матрица Cij на каждом процессорном эле менте.

В алгоритме Кэннона в оперативной памяти каждого процессор ного элемента также должны размещаться 3n2 / p чисел с плаваю щей точкой, однако объем пересылаемой между процессорными элементами информации намного меньше. Действительно, количест во чисел, пересылаемых за все время выполнения алгоритма (здесь также предполагается, что на каждом шаге алгоритма данные пере даются одновременно), равно s (этапов)*2(блока)* n2 / s 2 (размер одного блока)= 2n 2 / p.

шаг шаг шаг Алгоритм Фокса Пусть p s. Представим исходные матрицы A, B и результат C в виде наборов квадратных блоков размером n / s n / s.

A11 A12 A1s B11 B12... B1s C11 C12... C1s...

A21 A22... A2 s B21 B22... B2 s C21 C22... C2 s,..............................

......

As1 As 2... Ass Bs1 Bs 2... Bss C s1 C s 2... Css s Cij Aim Bmj, i, j 1,.., s.

m Для организации параллельных вычислений предположим, что процессоры образуют логическую прямоугольную решетку s s (pij, i, j 1,..., s).

В соответствии с алгоритмом Фокса каждый процессорный эле мент pij отвечает за вычисление Cij. В ходе вычислительного про цесса на каждом pij размещаются вычисляемое Cij, исходная матри ' ' ца Aij, Aij и Bij – блоки матриц A, B, получаемые в ходе выполнения алгоритма.

На этапе инициализации производится рассылка на каждый вы числительный узел pij блоков Aij, Bij и присваивание нулевых зна чений элементам матрицы Cij. При проведении вычислений на каж дой итерации l (1 l s ) выполняются следующие действия:

- для каждого ряда процессорной решетки i (1 i s ) определя ется значение индекса m по формуле m i l 1 mod s 1 ;

- блок Aim пересылается на все процессорные элементы ряда i и производится присваивание Aij Aim, j 1,..., s ;

' ' - полученные в результате пересылок матрицы Aij, Bij каждого процессорного элемента pij перемножаются и прибавляются к мат рице Cij :

' ' Cij Cij Aij Bij ;

' - матрицы Bij каждого процессорного элемента pij пересылаются pi 1 j, являющимися соседями сверху в колонках процессорной ре шетки (ПЭ p1 j отправляет данные ПЭ psj ).

шаг шаг шаг Приведенный параллельный метод матричного умножения в слу чае, если n нацело делится на s, приводит к равномерной загрузке процессоров:

t t n3 ;

n Tp tmult tadd s mult add s p tmult tadd n3 p.

Sp (2.17) tmult tadd n3 / p Однако блочное представление матриц в рассмотренных выше алгоритмах умножения матриц приводит к некоторому увеличению объёма пересылаемых между процессорными элементами данных.

Так, в алгоритме Фокса на каждой итерации этапа вычислений на pij процессорный элемент дополнительно передаётся два массива дан ных общим объёмом 2(n / s) 2. Тогда на каждый вычислительный узел необходимо отправить (без учета затрат на размещение блоков при подготовке к выполнению алгоритма) 2 n Vij 2( n / s) 2 s 2n2 / p p данных. Объем пересылаемых данных может быть уменьшен, если использовать строковое для A и столбцовое для B блочное распре деление данных по процессорным элементам.

Блочное представление матриц при параллельном вычислении матричных произведений имеет ряд преимуществ:

- пересылка данных является распределённой во времени, что по зволяет совмещать вычисления и передачу данных;

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

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

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

2.4. Библиотека базовых подпрограмм линейной алгебры Для реализации операций линейной алгебры на таких вычисли тельных системах, как CrayT3E, IBM SP2, Intel Paragon, SGI Power Challenge, TM CM-5, Linux-кластерах или любой другой системе, на которой установлена библиотека MPI (Message Passing Interface) или PVM (Parallel Virtual Machine), разработана специализированная библиотека PBLAS (основные параллельные подпрограммы линей ной алгебры – Parallel Basic Linear Algebra Subroutines).

Библиотека PBLAS включает библиотеку подпрограмм BLAS (Basic Linear Algebra Subroutines) для осуществления базовых опера ций линейной алгебры на серийных компьютерах и библиотеку BLACS (Basic Linear Algebra Communication Subroutines) для органи зации параллельных вычислений. Библиотека PBLAS является ос новной компонентой библиотеки ScaLAPACK (Scalable LAPACK), предназначенной для решения отдельной задачи линейной алгебры, например, решения системы линейных алгебраических уравнений, нахождения собственных значений вещественной симметричной матрицы, LU-разложения матрицы, приведения вещественной сим метричной матрицы к трехдиагональному виду [6,7].

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

Эти подпрограммы учитывают иерархию памяти современных про цессоров и делятся на три уровня. В библиотеке BLAS первого уров ня (BLAS1) собраны подпрограммы для операций над векторами или парами векторов. Они выполняют O (n) арифметических операций, результатом которых является либо вектор, либо число. Примерами подпрограмм BLAS1 являются: операция saxpy, скалярное произве дение, вычисление нормы вектора (табл. 2.1).

Таблица 2. Примеры подпрограмм библиотеки BLAS Уровень Кол-во Примеры Функция операций BLAS1 Saxpy Скалярвектор+вектор O (n ) Sdot Скалярное произведение Snrm2 Евклидовая норма вектора BLAS2 Sgemv Произведение матрицы на вектор O (n 2 ) Strsv Решение СЛАУ с треугольной матрицей BLAS3 Sgemm Произведение матриц O (n 3 ) Strsm Решение нескольких систем с треугольными матрицами BLAS второго уровня (BLAS2) объединяет подпрограммы для матрично-векторных операций. Эта часть библиотеки BLAS спроек тирована таким образом, чтобы оптимизировать использование дан ных кэш-памяти и уменьшить количество обменных операций между быстрой и медленной памятью. BLAS2 совершает операции главным образом с матрицей (двумерным массивом) и вектором (или векто рами), результатом которых являются матрица или вектор. Если мас сив имеет размерность n n, то необходимое количество операций оценивается как O (n 2 ).

BLAS третьего уровня (BLAS3) предназначена для реализации матричных операций и учитывает все преимущества современной архитектуры компьютера, в значительной степени минимизирует в вычислительном алгоритме отношение числа операций с плавающей точкой к количеству обращений к быстрой памяти. Данные, разме щенные в кэш-памяти или локальной памяти, оптимально использу ются за счет разбиения матриц на более мелкие блоки. Матричные операции осуществляются на основе блочного представления мат риц. BLAS3 оперирует с парами или тройками матриц. В результате получается матрица. Примерами подпрограмм из библиотеки BLAS являются подпрограмма умножения матриц и подпрограмма реше ния нескольких систем с треугольными матрицами (табл. 2.1).

Поскольку библиотека BLAS включает высокоэффективные под программы, может применяться в различных вычислительных сис темах и широко используется мировым научным сообществом, в на стоящее время она является основным средством при разработке вы сококачественного программного обеспечения для решения задач вычислительной линейной алгебры, например пакетов LINPACK и LAPACK. Библиотека BLAS распространяется свободно (http://www.netlib.org/blas).

BLACS – библиотека для обеспечения параллельных вычислений при решении задач линейной алгебры, использующая стандарт пере дачи сообщений. Эта библиотека организует параллельные вычисле ния на одно- или двумерной логической процессорной сетке, где ка ждый процесс обрабатывает части матриц или векторов. BLACS включает подпрограммы для синхронного получения или отправки сообщений, чтобы обеспечить передачу матриц или их блоков с од ного процесса другому, рассылку блоков матриц нескольким процес сам или глобальную редукцию данных (суммирование, определение максимального или минимального значения). В этой библиотеке также есть подпрограммы для построения и модификации логиче ской процессорной сетки. В настоящее время BLACS эффективно используется на многопроцессорных вычислительных системах, где установлены библиотеки передачи сообщений, например MPI или PVM.

3. ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ 3.1. Решение систем линейных уравнений с заполненными матрицами методом исключения Гаусса Рассмотрим систему линейных алгебраических уравнений Ax b (3.1) с невырожденной матрицей A размера n n. Будем считать матрицу A – плотной, т.е. содержащей относительно небольшое число нуле вых элементов.

Одним из прямых методов решения линейных систем является метод исключения Гаусса. Идея этого метода состоит в том, что мат рица A сначала упрощается – приводится эквивалентными преобра зованиями к треугольной или диагональной, а затем решается систе ма с упрощенной матрицей. Наиболее известной формой исключения Гаусса является та, в которой расширенная матрица системы линей ных уравнений приводится к верхнетреугольному виду путём вычи тания одних уравнений, умноженных на подходящие числа, из дру гих уравнений. Полученная треугольная система решается с помо щью обратной подстановки.

Применение метода исключения Гаусса основывается на LU теореме. Пусть дана квадратная матрица A порядка n, и пусть Ak – главный минор матрицы, составленный из первых k строк и столб цов. Если det Ak 0 для k 1, 2,..., n, то тогда существуют единст венная нижняя треугольная матрица L с единичными элементами на главной диагонали и единственная верхняя треугольная матрица U такие, что A LU.

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

LUx b. (3.2) Тогда исходная система (3.1) сводится к последовательному ре шению двух систем:

- решая систему с нижней треугольной матрицей Ly b прямой подстановкой, получим вектор y ;

- затем, решая систему с верхней треугольной матрицей Ux y обратной подстановкой, находим решение x исходной системы уравнений (3.1).

Рассмотрим вначале LU-разложение, поскольку проведение этой операции требует значительных временных затрат. Ниже представ лен псевдокод LU-факторизации, где значения диагональных и над диагональных коэффициентов матрицы U в ходе вычислительного процесса переписываются в соответствующие элементы матрицы A.

Пусть aii 0, i 1,..., n.

do k=1,n- do i=k+1,n lik = aik/akk end do do j=k+1,n do i=k+1,n aij = aij – lik*akj end do end do end do В цикле j производится вычитание k -й строки матрицы A, ум ноженной на соответствующее число, из расположенных ниже строк.

Правая часть b системы (3.1) также может быть добавлена к матри це A ( n +1 столбец) и обрабатываться в ходе приведения к треуголь ному виду. В общем случае может потребоваться перестановка строк для обеспечения корректности и вычислительной устойчивости LU разложения (выбор главного элемента в случае, когда aii 0 или aii 0 ).

В методе исключения Гаусса выполняется около n3 / 3 парных операций сложения и умножения, а также n 2 / 2 делений. Пренебре гая членом меньшего порядка, получим для последовательных вы числений T1 tmult t add n3 / 3.

LU-разложение (метод исключения Гаусса) представляется в виде тройных вложенных циклов, в которых рассчитанные элементы мат рицы U переприсваивают элементам исходной матрицы A :

do _ do _ do aij = aij – (aik/akk)/akj end do end do end do Индексы i, j, k в циклах могут быть взяты в любом порядке. Все го возможно 3! = 6 вариантов задания порядка следования циклов, как и для рассмотренного в предыдущей главе в п.2.3 умножения квадратных матриц. Шесть различных форм следования ijk -циклов отличаются способом доступа к элементам матриц A и L lik, ко торый может существенно повлиять на производительность алгорит ма в зависимости от архитектуры компьютера. Формы jki и kji столбцово-ориентированы: они обращаются к столбцам матриц A и L. В формах kij и ikj происходят обращения к строкам матрицы A и лишь к отдельным элементам L. Для оставшихся форм ( ijk и jik ) характерен смешанный доступ – к столбцам в A и к строкам в L.

По-видимому, наиболее оптимальными для параллельных реализа ций являются формы kij и kji, которые отличаются только в форме доступа к элементам матрицы A – по строкам или по столбцам соот ветственно.

Применим общий подход построения параллельных процедур ре шения задач для алгоритма LU-разложения.

1. Декомпозиция. Выберем в качестве отдельных блоков (фунда ментальных подзадач) алгоритма метода исключения Гаусса сле дующие: для i, j 1,..., n каждая фундаментальная задача ( i, j ) по значению aij вычисляет и запоминает значения uij, i j.

lij, i j Все такие фундаментальные задачи образуют двумерный массив n2 подзадач (рис. 3.1). Нулевые элементы матриц U [uij ] и L, а также единичные элементы на главной диагонали матрицы L не рас считываются.

Рис. 3.1. Фундаментальные подзадачи и их коммуникационная схема (n = 8) 2. Организация коммуникаций. Необходимые коммуникации для связи фундаментальных подзадач для обеспечения всего цикла вы числений представляет следующий псевдокод для задачи ( i, j ):

do k=1,min(i,j)- получить akj получить lik aij = aij – lik*akj end do if ij then разослать aij подзадачам (k,j), k=i+1,…,n else получить ajj lij = aij/ajj разослать lij подзадачам (i,k), k=j+1,…,n end if Из рис. 3.1 видно, что каждая фундаментальная подзадача ( i, j ) связана с подзадачами (1, j ),...,( n, j ) и (i,1),...,(i, n).

3. Укрупнение. Для имеющего n n массива фундаментальных подзадач могут быть использованы следующие стратегии получения p ( p n 2 ) составных подзадач:

- одномерное укрупнение: объединяются n/p строк (или столб цов) двумерного массива фундаментальных подзадач;

- двумерное укрупнение: одна составная подзадача включает в себя (n / p ) ( n / p ) фундаментальных подзадач.

4. Распределение составных подзадач по процессорным элемен там производится с учетом архитектуры многопроцессорной вычис лительной системы.

Сначала рассмотрим одномерное укрупнение, когда составная подзадача объединяет n/p строк массива фундаментальных подзадач (рис. 3.2).

Рис.3.2. Одномерное укрупнение строк массива фундаментальных подзадач В этом случае нет необходимости рассылать множители lij (коэф фициенты матрицы L ) по строкам, так как любая строка заданной матрицы целиком принадлежит одной составной подзадаче. Также следует заметить, что при данном способе укрупнения (когда строки последовательно по своему номеру распределяются составным под задачам) нет никакого параллелизма при расчете новых значений любой заданной строки. После пересчета элементов в строке необхо дима рассылка полученных значений, чтобы обеспечить вычисления для каждой строки нижней составной подзадачи (рис. 3.2). Рассмот ренная параллельная процедура LU-разложения может быть пред ставлена следующим псевдокодом:

do k = 1,n- if k myrows then переслать {akj: kjn} другим составным подзадачам else получить {akj: kjn} end if do i myrows, ik lik = aik/akk end do do j = k+1,n do i myrows, ik aij = aij – lik*akj end do end do end do Здесь myrows – множество индексных значений номеров строк, рас пределенных одной составной подзадаче.

В этом алгоритме на каждом шаге k каждой задаче требуется около (n k ) 2 / p арифметических операций. Поэтому без учета коммуникационных затрат получим n Tpcomp tmult tadd n k / p tmult tadd n 3 / 3 p.

k Количество данных, которые необходимо переслать на каждом шаге, около ( n k ). Тогда n Tpcomm tcomm n k tcomm n 2 / 2.

k Здесь tcomm – время на передачу одного числа. Общее время для вы числения LU-разложения на p-процессорной вычислительной систе ме будет Tp Tpcomp Tpcomm tmult t add n3 / 3 p tcomm n 2 / и tmult tadd n3 / T1 p. (3.3) Sp Tp tmult tadd n3 / 3 p tcomm n 2 / 2 3 ptcomm 2n tmult tadd Чтобы уяснить особенности построенной процедуры для парал лельного LU-разложения, рассмотрим случай, когда p n. Тогда i -я строка матрицы A хранится в i -м процессорном элементе. На пер вом шаге вычислительного процесса первая строка рассылается на все вычислительные узлы, а процессорные элементы с номерами 2,3,..., n выполняют параллельные вычисления:

li1 ai1 / a11, aij li1a1 j, j 2,..., n.

aij На втором шаге процессорный элемент с номером 2 рассылает вто рую строку преобразованной матрицы A процессорным элементам 3, 4,..., n, на которых вновь проводятся параллельные вычисления:

li 2 ai2 / a22, aij aij li 2 a2 j, j 3,..., n.

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

Таким образом, для одномерного строчно-ориентированного ук рупнения в алгоритме LU-разложения характерны следующие осо бенности:

1. Каждая составная подзадача считается выполненной, как толь ко заканчивается обработка последней строки распределенной на процессорный элемент части матрицы A. Если составная подзадача объединяет строки матрицы A с последовательной нумерацией, то обрабатывающий строку процессорный элемент завершает свою ра боту задолго до окончания общего вычислительного процесса.

2. Обработка строк в порядке увеличения их номера i требует существенно меньшего объема вычислительной работы (по сравне нию, например, с обработкой строки с номером 1), что ведет к нерав номерной загрузке процессоров.

3. Балансировка загрузки процессоров может быть улучшена, если строки матрицы A для составных подзадач объединять циклически, т.е. строка i назначается процессорному элементу с номером i mod p, i 1,..., n (рис. 3.3).

строки: строки: строки: строки:

1 2 m p p+1 p+2 p+m 2p … … … … … … (k-1)p+1 (k-1)p+2 (k-1)p+m kp ПЭ 1 ПЭ 2 ПЭ m ПЭ p Рис. 3.3. Циклическая строчная схема распределения фундаментальных подзадач по процессорным элементам (k = n/p) Проблема эффективности загрузки процессоров для этой схемы несколько теряет остроту. Так, процессорный элемент с номером (ПЭ 1) производит вычисления при p n практически до оконча ния процесса факторизации матрицы A. Тем не менее некоторая не равномерность загрузки процессоров при такой схеме сохраняется:

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

Также можно использовать другие способы распределения дан ных по процессорным элементам. Например, циклическое размеще ние блоков матрицы A или применение слоистой схемы с отраже ниями (рис. 3.4). Для этой схемы первые p строк матрицы A рас пределяются между p процессорными элементами в естественном порядке, следующие p строк размещаются в обратном порядке и т.д.

строки:

строки: строки: строки:

p 1 2 m p+ 2p 2p-1 2p-m … … … … … … kp (k-1)p+1 (k-1)p+2 (k-1)p+m ПЭ p ПЭ 1 ПЭ 2 ПЭ m Рис. 3.4. Слоистая строчная схема с отражениями, k = n/p – нечетное 4. Производительность параллельных вычислений может быть увеличена за счет совмещения коммуникационных операций и вы числений. Рассмотрим, как это можно выполнить для циклической схемы распределения данных (рис. 3.3). Пусть процесс обработки строк ведется в порядке увеличения их номера. Тогда на первом шаге каждый процессорный элемент, получив от ПЭ 1 первую строку, сначала вычисляет множители li1, а затем модифицирует соответст вующие строки матрицы A. На втором шаге процессорный элемент, хранящий вторую строку преобразованной матрицы, должен пере слать ее на остальные вычислительные узлы. Здесь возникает за держка вычислительного процесса. Выходом из этой ситуации явля ется немедленная рассылка обработанной второй строки, пока все процессоры будут заняты модификациями первого шага. Подобным образом можно организовать параллельные вычисления и на других шагах: на m -ом шаге ( m +1)-я строка рассылается сразу после того, как закончится ее обработка. Такая стратегия осуществления парал лельных вычислений называется опережающей рассылкой. Ее эф фективность в значительной степени определяется топологией меж процессорных соединений МВС с распределенной памятью. Страте гию опережающей рассылки можно использовать и для систем с раз деляемой памятью. Однако в силу особенностей архитектуры МВС в данном случае эта стратегия носит название опережающего вычис ления, поскольку меняется порядок действий. Как только закончится модификация ( m +1)-й строки, она маркируется признаком «готова».

Тогда другие процессоры, завершив свою работу на m -ом шаге, мо гут сразу приступать к следующему ( m +1)-му шагу. Маркировка дает возможность осуществлять необходимую синхронизацию рабо ты процессоров без неоправданных задержек. Часто опережающую рассылку и опережающее вычисление называют конвейеризацией.

Рис. 3.5. Столбцово-ориентированное укрупнение фундаментальных подзадач для LU-разложения Альтернативой строкового укрупнения фундаментальных подза дач является их объединение по столбцам (рис. 3.5). Заметим, что в этом случае нет необходимости в пересылке строк матрицы A, по скольку все необходимые для вычислений данные принадлежат од ной составной подзадаче или процессорному элементу. Однако при таком способе распределения данных нет явного параллелизма при вычислении множителей lij и модификации элементов в столбцах матрицы. Также очевидна необходимость организации рассылки множителей lij подзадачам, расположенным справа (рис. 3.5), для модификации элементов матрицы A.

Запишем псевдокод параллельной процедуры LU-разложения при распределении матрицы A по столбцам в процессорных элементах.

do k = 1,n- if k mycols then do i = k+1,n lik = aik/akk end do переслать {lik: kin} другим составным подзадачам else получить {lik: kin} end if do j mycols, jk do i = k+1,n aij = aij – lik*akj end do end do end do Здесь mycols – множество индексных значений номеров столбцов, распределенных одной составной подзадаче (процессорному элемен ту). Параметры ускорения и эффективности данного алгоритма весь ма похожи на соответствующие выражения для строчно ориентированного распределения данных (3.3).

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

1. Каждая составная подзадача (укрупненный блок) считается вы полненной, как только закончится обработка ее последнего столбца.

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

2. Расчет множителей и модификация столбцовых элементов мат рицы требуют значительно меньшего объема вычислений при увели чении номера столбца. Для повышения согласованности работы про цессоров и балансировки загрузки можно использовать циклическое распределение столбцов матрицы A по процессорным элементам:

например, столбцы матрицы j 1,..., n размещать на процессорные элементы, номер которых определяется как j mod p. Кроме того, мо гут применяться другие способы распределения данных – блочно циклический или слоистая столбцовая схема с отражениями.

3. Производительность построенной параллельной процедуры может быть улучшена за счет совмещения обменных операций с вы числениями. На m -ом шаге каждый процессор выполняет модифи кацию своей части необработанной матрицы, прежде чем осущест вится переход к следующему ( m +1)-му шагу алгоритма. Однако процессорный элемент, на котором размещается m +1 столбец мат рицы, может провести вычисления множителей для шага m +1 и их рассылку на другие ПЭ перед завершением всех модификаций эле ментов матрицы на m -ом шаге. Стратегия опережающей рассылки позволяет процессорам начинать на каждом шаге вычислительную работу раньше по сравнению с обычным подходом, когда они ожи дали получение необходимых для вычисления данных от других ПЭ.

Другим способом укрупнения является объединение n / p n / p фундаментальных подзадач в один блок (рис. 3.6).

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

Итоговый алгоритм является блочной версией мелкозернистого ал горитма.

Параллельный псевдокод алгоритма LU-разложения для двумер ного укрупнения можно представить следующим образом:

do k = 1,n- if k myrows then разослать {akj: j mycols, jk} другим подзадачам в столбце этой подзадачи else получить {akj: j mycols, jk} end if if k mycols then do i myrows, ik lik = aik/akk end do разослать {lik: i myrows, ik} другим составным подзадачам в строке этой подзадачи else получить {lik: i myrows, ik} end if do j mycols, jk do i = myrows, ik aij = aij – lik*akj end do end do end do Рис. 3.6. Двумерное укрупнение мелкозернистых подзадач На каждом шаге k количество операций, необходимых для мо дификации элементов матрицы каждой составной подзадачи, опре деляется как (n k ) 2 / p. Тогда общее время на проведение парал лельных вычислений без учета обменных операций составит n Tcomp tmult tadd n k / p t mult t add n3 / p.

k Количество рассылок данных на шаге k вдоль каждой строки или столбца составных подзадач (рис. 3.6) около (n k ) / p. Тогда для двумерной сетки общее время на обменные операции можно запи n p tcomm n 2 / сать как Tcomm 2 tcomm n k / p.

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

Общее время оценивается Tp tmult tadd n3 / (3 p ) tcomm n 2 / p и tmult tadd n3 / Sp tmult tadd n3 / (3 p) tcomm n2 / p (3.4) p.

3tcomm p tmult tadd n Перечислим некоторые особенности параллельной реализации LU-разложения при двумерном укрупнении мелкозернистых подза дач:

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

2. При вычислении множителей и модификации элементов матриц все меньшее количество арифметических операций необходимо при увеличении номеров строк и столбцов для последовательного рас пределения элементов матриц по подзадачам. Согласованность рабо ты процессоров и балансировка загрузки процессорных элементов может быть улучшена при использовании циклического способа рас пределения данных, когда элемент матрицы aij, i, j 1,..., n присваи вается составной подзадаче с номером ( i mod p, j mod p, рис. 3.6).

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

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

Частичный выбор главного элемента Чтобы обеспечить корректность и устойчивость LU-разложения, часто применяется схема частичного выбора главного элемента, в результате которой в непреобразованной части матрицы A произво дится перестановка строк таким образом, чтобы в ведущем столбце находился наибольший по модулю элемент матрицы. При такой стратегии множители lik по абсолютной величине не будут превы шать единицы, что способствует уменьшению влияния ошибки ок ругления на конечный результат.

Использование схемы частичного выбора главного элемента для обеспечения устойчивости процесса усложняет параллельную реали зацию исключения Гаусса и оказывает значительное влияние на про изводительность параллельных вычислений. При одномерной столб цовой декомпозиции поиск главного элемента не требует обмена данными между процессорными элементами – он ведется каким-то одним процессором. Однако в этом случае остальные процессоры простаивают. Выход из такой ситуации может быть обеспечен путем опережающего вычисления и рассылки. На k -м шаге поиск главного элемента в k +1-м столбце сразу же начинается после модификации этого столбца. Как только будет найдена ведущая строка, эта инфор мация сразу же передается другим процессорам, которые могут па раллельно выполнить перестановку строк.

В заключение заметим, что систему линейных алгебраических уравнений Ax b можно решить, пользуясь формулами Краме ра ( det A 0 ) xi i /, i 1, 2,..., n, (3.5) где i – определитель, получающийся из матрицы A заменой i-го столбца столбцом правых частей системы. Непосредственное вычис ление определителя матрицы A по формулам линейной алгебры требует O(n!) операций, что при n 1 не экономично. Поэтому для вычисления можно воспользоваться описанным выше алгоритмом прямого хода метода Гаусса (не выполняя действий со столбцом правой части системы), в результате которого будет получена верх нетреугольная матрица U [ui, j ]1 за O(n3) операций, определитель n которой по LU-теореме равен определителю матрицы. Все опре делители из (3.5) можно вычислить за O(n4) операций. Отметим, что метод Гаусса можно применить и для нахождения обратной матри цы, учитывая, что A1 A E.

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

3.2. Решение систем с треугольными матрицами После выполнения LU-разложения (3.2) необходимо решать сис темы линейных уравнений с треугольными матрицами: Ly b и Ux y. Заметим, что большинство прямых методов решения линей ных систем с заполненными матрицами приводят матрицу A к тре угольному виду, а затем ищется решение полученной эквивалентной системы с треугольной матрицей. Кроме того, системы с треуголь ными матрицами часто используются для предобусловливания в ите рационных методах решения СЛАУ.

Для системы с нижней треугольной матрицей Lx b решение может быть получено в результате прямой подстановки:

i xi bi lij x j / lii, i 1,..., n. (3.6) j Псевдокод этого алгоритма можно записать следующим образом:

do i = 1,n xi = bi/lii do j = i+1,n bj = bj – lji * xi end do end do При решении системы с верхней треугольной матрицей Ux b требуется применять обратную подстановку:

n xi bi uij x j / uii, i n,...,1. (3.7) j i В этом случае псевдокод имеет вид do i = n, xi = bi/uii do j = 1,i- bj = bj – uji * xi end do end do Заметим, что представленный выше псевдокод для решения сис темы Ux b полезен, если матрица U хранится в памяти компьюте ра по столбцам. Этот алгоритм называется столбцовым алгоритмом (или алгоритмом векторных сумм). Альтернативный алгоритм, ко торый носит название алгоритма скалярных произведений, имеет преимущество, когда матрица U хранится по строкам. Он задает ска лярные произведения с длинами векторов, меняющимися от 1 до n 1, и n делений чисел.

do i = n, do j = i+1,n bi = bi – uij * xj end do xi = bi/uii end do Анализируя эти два алгоритма, можно отметить, что в столбцовом алгоритме обратной подстановки производится немедленное (а не отложенное) определение значения i -й неизвестной. Однако какой из этих алгоритмов следует выбрать, диктуется, главным образом, способом хранения матрицы U, если он был определен в процессе LU-разложения.

Подсчитаем количество арифметических операций сложения или умножения для получения решения по (3.6) или (3.7):

1+2+…+ n –1= n( n 1) / 2 n2 / 2.

Пренебрегая временем выполнения n операций деления в (3.6) или (3.7), можно записать T1 tadd tmult n 2 / 2. (3.8) Далее рассмотрим только нижнетреугольную систему Lx b, так как решение верхнетреугольной системы организуется аналогичным образом.

Разберем, как строятся параллельные алгоритмы для решения систем с нижнетреугольными матрицами.

Рис. 3.7. Декомпозиция и проектирование коммуникаций для алгоритма решения системы с треугольной матрицей На этапе декомпозиции в качестве фундаментальных подзадач будем выбирать следующие (см. рис. 3.7):

- для i 2,..., n;

j 1,..., i 1 ( i, j )-я подзадача запоминает lij, по лучает x j и вычисляет произведение lij x j ;

- для i 1,..., n ( i, i )-я подзадача хранит значения lii и bi, собирает i произведения lij x j, вычисляет сумму ti lij x j и вычисляет и запо j минает xi (bi ti ) / lii.

Таким образом, декомпозиция дает двумерный треугольный мас сив n( n 1) / 2 мелкозернистых подзадач.

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

- для j 1,..., n 1 ( j, j )-я подзадача одновременно передает зна чение x j подзадачам ( i, j ), i j 1,..., n ;

- далее производится пересылка произведений lij x j от подзадач ( i, j ), j 1,..., i 1, подзадаче ( i, i ).

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

Для ( i, j )-й подзадачи ( i j ) if i=j then if i1 then получить{t} от подзадач (i,j), j=1,…,i-1 (t-сумма ti) else t= end xi=(bi-t)/lii разослать {xi} подзадачам (k,i), k=i+1,…,n else получить {xj} t=l(i,j)*x(j) переслать {t} для суммирования подзадаче (i,i) end if Согласно этому фрагменту параллельной программы сначала должно быть найдено значение неизвестного x1, т.е. решена подза дача (1,1). На следующем шаге это значение пересылается подзада чам (2,1) и (2,2);

одновременно вычисляются произведения;

полу ченные каждой подзадачей значения t затем пересылаются подзада че (2,2), где рассчитывается значение x2. Далее процесс продолжает ся, вычисляется следующее неизвестное, и видно, что с увеличением индекса неизвестных все больше вычислений может проводиться одновременно, однако также увеличивается и количество коммуни каций между подзадачами. Чтобы уменьшить затраты на передачу значений между подзадачами, необходимо провести объединение мелкозернистых подзадач.

Укрупнение для имеющегося двумерного массива фундаменталь ных мелкозернистых подзадач может быть выполнено следующими способами:

- одномерным укрупнением n / p строк (или столбцов);

- двумерным укрупнением n / p n / p подзадач.

В итоге получается p крупнозернистых подзадач.

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

Рассмотрим одномерное укрупнение по строкам мелкозернистых подзадач (рис. 3.8).

Рис. 3.8. Одномерное строковое укрупнение мелкозернистых подзадач Видно, что такой способ укрупнения не требует передачи данных для вычисления суммы значений t по строкам, поскольку строка матрицы мелкозернистых подзадач целиком принадлежит укрупнен ной подзадаче. Заметим, что нет никакого параллелизма при вычис лении этих сумм. Кроме того, передача вычисленных значений неиз вестных другим укрупненным подзадачам все еще необходима.

Псевдокод параллельной программы можно записать следующим образом:

do j=1,n if j myrows then xj=bj/ljj разослать {xj} остальным подзадачам else получить {xj} end if do i=myrows, ij bi = bi – lij*xj end do end do (а) (б) Рис. 3.9. Циклическая схема (а) и схема с отражениями (б) распределения строк мелкозернистых подзадач В соответствии с принятой схемой укрупнения можно отметить следующие особенности полученной параллельной программы:

1. Выполнение каждой укрупненной подзадачи заканчивается по сле решения последнего уравнения, принадлежащего этой подзадаче.

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

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

4. Согласованность работы процессоров и балансировка их за грузки может быть улучшена за счет объединения строк в блоки в циклическом порядке, когда i -я строка фундаментальных подзадач назначается укрупненной подзадаче с номером i mod p (рис. 3.9,а).

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

3.9,б).

Подсчитаем временные затраты параллельного алгоритма. Время, необходимое для проведения вычислений по рассмотренному алго ритму, оценим как n Tpcomp n tdiv tadd tmult, 2p здесь tdiv – время операции деления двух чисел с плавающей точкой, а время, затрачиваемое на рассылку рассчитанных значений неиз вестных, представим в виде n p Tpcomm tcomm.

p В итоге n p n Tp tadd tmult tcomm. (3.9) 2p p Проведем оценку ускорения построенного алгоритма по форму лам (3.8) и (3.9) в предположении, что обменные операции между узлами перекрываются, т.е. происходят одновременно:

T p Sp 1. (3.10) 2tcomm p Tp tadd tmult n Из (3.10) следует, что повышение эффективности построенной параллельной программы может быть достигнуто при уменьшении доли временных затрат на обеспечение пересылки данных между вычислительными узлами. Одним из путей решения этой проблемы является использование стратегии опережающей рассылки, т.е. когда рассылка данных совмещается с вычислениями. Для рассматривае мой задачи на j -м шаге (см. представленный выше алгоритм) круп нозернистая подзадача, включающая строку (уравнение) j 1, может рассчитать значение x j 1 и разослать его другим подзадачам перед осуществлением окончательных вычислений с использованием x j.

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

Например, при циклической схеме распределения (рис. 3.9,а) на МВС с топологией «кольцо» имеет место практически полное пере крытие обменов, в то время как топология «гиперкуб» позволяет де лать существенно меньшее количество перекрытий.

Рассмотрим одномерное укрупнение фундаментальных мелкозер нистых подзадач по столбцам. Для такого объединения подзадач в p блоков (рис. 3.10) характерными являются следующие особенности:

1. Нет необходимости рассылать вычисленные значения неиз вестных x j по вертикали, поскольку каждый необходимый для расчетов столбец матрицы целиком принадлежит только одной ук рупненной подзадаче.

2. Нет никакого параллелизма при расчете произведений lij x j.

l x 3. Для вычисления суммы произведений потребуются пе ij j j ресылки между укрупненными блоками по горизонтали.

4. Каждая укрупненная подзадача начинает выполняться только тогда, когда начнет вычисляться неизвестная x j ее первого столбца.

5. Если каждая крупнозернистая подзадача объединяет последова тельные блоки столбцов, она может оставаться неактивной большую часть процесса решения задачи.

6. Количество произведений, осуществляемых с одной компонен той вектора x, уменьшается с увеличением номера столбца.

7. Равномерная загрузка процессорных элементов может быть улучшена за счет назначения столбцов блокам по циклической схе ме, когда столбец j распределяется укрупненному блоку с номером j mod p. Кроме того, могут использоваться другие схемы, например блочно-циклическая или схема с отражениями.

Рис. 3.10. Столбцово-ориентированная схема распределения мелкозернистых подзадач по блокам Оценки показывают, что для одномерного столбцового укрупне ния время выполнения параллельного алгоритма, в сущности, точно такое же, как и для одномерного строкового и аналогично зависит от перекрытия последовательных шагов алгоритма. Ускорение парал лельного алгоритма может быть улучшено за счет применения стра тегии опережающих вычислений.

do i=1,n t= do jmycols, ji, t=t+lijxj end do if imycols then получить {t} и провести суммирование полученных значений xi=(bi-t)/lii else отправить {t} для вычисления суммы end if end do Параллелизм рассмотренных выше алгоритмов получается из внутреннего цикла, выполнение которого распределяется по укруп ненным задачам, при этом внешний цикл выполняется последова тельно. Отметим, что эти параллельные алгоритмы вычисляют толь ко по одной компоненте решения x. Для повышения эффективности можно применить конвейерный способ, в котором используется яв ный параллелизм и для внешнего цикла, что позволяет рассчитывать несколько компонент вектора решения одновременно. Такие алго ритмы получили название волновых.

Действительно, рассмотрим одномерную схему укрупнения по столбцам (рис. 3.11). Кажется, что этот алгоритм не имеет никакого параллелизма: после того, как укрупненная подзадача, располагаю щая столбцом j, вычислит x j, окончательное обновление правой части не может быть разделено с другими подзадачами для паралле лизации, так как они не имеют доступа к столбцу j. Вместо того, что бы выполнять все такие обновления немедленно, подзадача, распола гающая столбцом j, могла бы рассчитать только, скажем, s компо нентов обновляемого вектора и пересылать (направлять) их подзада че, имеющей столбец j 1, перед продолжением работы со следую щими s компонентами обновляемого вектора. По мере получения первых s компонентов обновляемого вектора подзадача, содержа щая столбец j 1, может вычислять x j 1, проводить дальнейшие об новления по собственному вхождению в другие подзадачи и т.д.

Рис. 3.11. Одномерный столбцовый волновой алгоритм Чтобы формализовать волновой столбцовый алгоритм, введем вектор z, в котором будут накапливаться обновления вектора пра вых частей b, и определим сегмент, в котором содержится, по край ней мере, s последовательных компонентов z.

do jmycols do k=1, #segments получить {segment} if k=1 then xj=(bj-zj)/ljj segment=segment-{zj} end if do zisegment zi=zi+lijxj end do if |segment|0 then отправить {segment} подзадаче со столбцом j+ end if end do end do В зависимости от размера сегмента, расположения столбца, отно шения скорости обменов к скорости вычислений и т.п. все подзадачи возможно загрузить одновременной работой, если ее выполнять для различных компонентов решения. В приведенном выше псевдокоде Segment – это регулируемый параметр, который контролирует соот ношение временных затрат между обменами и параллельными вы числениями. «Первый» сегмент для данного столбца сдвигается на один элемент после того, как каждая компонента решения найдена, исчезает после s шагов, затем следующий сегмент становится «пер вым» сегментом и т.д. К концу выполнения алгоритма остается толь ко один сегмент, и он содержит только один элемент. Таким образом, в этом алгоритме снижается общее количество обменов.

Время, затрачиваемое на выполнение волнового столбцового ал горитма, может быть оценено следующим образом:

Tp tcomm tadd tmult n 2 np s s 1 p 2 / 2 p. (3.11) По мере того как длина сегмента s увеличивается, затраты на за пуск коммуникационных операций уменьшаются, но растут затраты на вычисления. Обратная картина имеет место для (3.11) при умень шении длины сегмента. Оптимальное значение длины сегмента s может быть определено из модели.

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

3.3. Решение линейных систем с трёхдиагональными матрицами Рассмотрим систему линейных алгебраических уравнений с трех диагональной матрицей следующего вида:

b0 c0 0 x0 f 0 0 a1 b1 c1 0 x1 f 0 0 a2 0 x2 f b2 c2. (3.12) 0 0 xn2 f n an2 bn2 cn an1 bn1 cn1 xn1 f n 0 0 0 bn xn f n an 0 0 Системы (3.12) образуют очень важный класс линейных алгеб раических уравнений. Они часто получаются в результате разност ных аппроксимаций дифференциальных краевых задач, а также при построении кубических сплайнов. Экономичными прямыми метода ми решения таких систем уравнений на компьютерах с последова тельной архитектурой являются специальные варианты метода ис ключения Гаусса – метод прогонки и методы редукции.

К сожалению, многие методы, которые оказались наиболее эф фективными в последовательных алгоритмах, не подходят для пря мого использования на параллельных компьютерах. Примером тако го метода является метод прогонки.

Метод прогонки Рассмотрим формулы метода прогонки для последовательных вы числений. Пусть необходимо решить следующую систему линейных уравнений:

ai xi 1 bi xi ci xi 1 fi, i 0,..., n;

(3.13) a0 0, cn 0.

Будем рассматривать только системы вида (3.13), которые имеют диагональное преобладание ( bi ai ci, i 0,..., n, причем хотя бы в одном из перечисленных неравенств обязательно должно выпол няться строгое неравенство).

Для решения системы (3.13) методом прогонки необходимо вы полнить:

- вычисление прогоночных коэффициентов («прямой ход» мето да) c f P0 0 ;

Q0 0 ;

b0 b ci f ai Qi, Qi i, i 1,..., n 1 ;

(3.14) Pi bi ai Pi 1 bi ai Pi ( Pi, Qi – прогоночные коэффициенты) - определение значений неизвестных («обратный ход» метода) f an Qn xn n ;

xi Pi xi 1 Qi, i n 1,..., 0. (3.15) bn an Pn Количество арифметических операций, необходимое для прове дения расчетов по формулам (3.14), (3.15), определяется как 6 n 1 2n 7 8n. То есть при увеличении размера системы (3.12) количество последовательных вычислений или время счета по по следовательной программе, реализующей метод прогонки, возрастет линейно.

Заметим, что формулы (3.14), (3.15) представляют не единствен ный способ применения метода прогонки. Например, если рассмат ривать вместо (3.15) рекуррентную формулу, в которой неизвестные линейной системы определяются в порядке увеличения индекса i, то получается:

- «прямой ход» – вычисление прогоночных коэффициентов Ri,Ti :

an f, Tn n ;

Rn bn bn ai f cT, Ti i i i 1 ;

i n 1,...,1 ;

(3.16) Ri bi ci Ri 1 bi ci Ri - «обратный ход» – вычисление неизвестных f c T x0 0 0 1 ;

xi Ri xi 1 Ti, i 1,.., n. (3.17) b0 c0 R Из формул (3.16) и (3.17) видно, что построенный вариант метода прогонки имеет такое же количество арифметических операций, что и (3.14), (3.15), а отличается только порядком изменения индексной переменной при вычислении прогоночных коэффициентов и неиз вестных. Оба рассмотренных метода решения системы линейных уравнений (3.12) являются последовательными, поскольку в них рас сматриваются рекуррентные формулы – как линейные (о способах распараллеливании которых указывается в главе 1), так и нелиней ные. При выборе варианта метода прогонки преимущество имеет тот, в котором первые прогоночные коэффициенты ( Pi или Ri ) име ют наименьшие по абсолютной величине значения, что в результате снизит влияние ошибки округления на конечный результат.

Если применять формулы (3.14), (3.15) для вычисления прогоноч ных коэффициентов Pi, Qi и неизвестных xi для значений ин Ri,Ti декса i 0,..., n / 2 1, а (3.16), (3.17) – для вычисления и xi при i n,..., n / 2 1, то можно получить вариант метода прогон ки, в котором вычисления на этапах прямого и обратного хода могут проводиться независимо и допускать распараллеливание построен ного алгоритма.


Рис. 3.12. Метод встречных прогонок, n = Действительно (см. рис. 3.12), для своих подобластей изменения индексной переменной i одновременно могут находиться прогоноч ные коэффициенты Pi, Qi и Ri,Ti. И если известно значение xn /2, то можно по формулам (3.15) и (3.17) одновременно и незави симо вычислить все неизвестные. Для равномерной загрузки процес соров необходимо, чтобы n было четным.

Такой вариант метода прогонки получил название «метод встреч ных прогонок». Заметим, что значение неизвестной xn /2 может быть найдено после первого этапа метода встречных прогонок – опреде ления прогоночных коэффициентов. Для этого в уравнении an /2 xn /2 1 bn /2 xn /2 cn /2 xn /2 1 f n / необходимо, пользуясь формулами (3.15) и (3.17), исключить xn /2 1.

Тогда получим f an /2 Qn /21 cn /2Tn /2 xn /2 n /2. (3.18) bn /2 an /2 Pn /2 1 cn / 2 Rn /2 Идея метода встречных прогонок была обобщена в работе Н.Н. Яненко с сотрудниками на случай многопроцессорной системы и получила название параметрической прогонки. Основу предлагае мого метода составляют:

- выделение части искомых неизвестных xi в качестве парамет рических (обозначим их x );

- построение системы линейных уравнений относительно этих не известных;

- решение этой системы и определение x ;

- нахождение по x остальных искомых неизвестных xi.

Для осуществления этой идеи система уравнений (3.12) разбива ется на подсистемы, количество которых будет совпадать с количе ством используемых процессорных элементов p. Причем декомпо зиция по данным осуществляется таким образом, что каждый -й процессорный элемент обрабатывает m 1 n / p 1 уравнений (см.

рис. 3.13), а для неизвестных системы вводятся локальные индексы для получения однородного алгоритма:

xi x j m x j, i 0,..., n;

j 0,..., m;

0,..., p 1.

Рис. 3.13. Пример декомпозиции неизвестных между тремя процессорами Кроме того, принимается, что xm1) x0 ) x, 1,..., p 1;

x0 x0 ;

x p xn.

( ( (3.19) Решение системы (3.12) будем искать в виде рекуррентных формул следующего вида:

x (j ) v (j ) x z (j ) x1 w(j ) ;

j 0,..., m ;

0,..., p 1. (3.20) v, z,w ( ) ( ) ( ) Здесь – решения следующих систем линейных j j j уравнений с трехдиагональной матрицей:

a j m v (j1 b j m v (j ) c j m v (j1 0;

j 1,..., m 1;

) ) (3.21) v0 ) 1, vm ) 0.

( ( a j m z (j 1 b j m z (j ) c j m z (j 1 0;

j 1,..., m 1;

) ) (3.22) z0 ) 0, zm ) 1.

( ( a j m w(j1 b j m w(j ) c j m w(j1 f j m ;

j 1,..., m 1;

) ) (3.23) w0 ) 0, wm ) 0.

( ( Системы (3.21) – (3.23) легко получаются из подстановки (3.20) в (3.13), решаются независимо для каждой -й подобласти обычным методом прогонки и допускают сокращение количества арифметиче ских операций за счет использования одной и той же трехдиагональ ной матрицы коэффициентов в (3.21) – (3.23).

После вычисления необходимых коэффициентов производится построение системы линейных уравнений для параметрических не известных x. Для этого рассматриваются неиспользованные ранее уравнения системы (3.13) при i 0, m,2m,..., pm, в которых с помо щью формулы (3.20) исключаются все «непараметрические» неиз вестные. В результате получим систему с трехдиагональной матри цей ( a0 0, cn 0 ):

am vm)1 x1 bm am zm)1 cm v1( ) x cm z1( ) x ( ( (3.24) f m am wm1) cm w1( ) ;

0,..., p.

( Решая эту систему обычным методом прогонки, найдем x. Да лее необходимо лишь воспользоваться формулой (3.20), чтобы опре делить остальные неизвестные задачи.

Хотя данный вычислительный алгоритм параллельного решения системы линейных алгебраических уравнений с трехдиагональной матрицей обладает высокой степенью параллелизма (вычисления по формулам (3.20) – (3.23) производятся с наивысшей для данной зада чи степенью параллелизма p ), тем не менее решение системы для параметрических неизвестных (3.24) выполняется последовательно только одним процессором. Кроме того, при использовании много процессорных вычислительных систем с распределенной памятью вычисление коэффициентов (3.24) и значений неизвестных по (3.20) также будет связано с дополнительными коммуникационными затра тами, обусловленными необходимой передачей данных из локальной памяти других процессоров.

Оценим временные затраты на получение решения по данному алгоритму. На первом этапе (решение систем (3.21) – (3.23)) каждому процессору требуется около 13 m арифметических операций, кото рые выполняются одновременно и независимо от других процессо ров. Далее решается система (3.24) одним из процессоров. Для рас чета коэффициентов системы (3.24) и нахождения ее решения требу ется около 18 p арифметических операций. Однако, если расчет ве дется на МВС с распределенной памятью, то потребуется доставить в локальную память этого процессора 6 p чисел и возвратить осталь ным процессорам полученное решение системы (3.24) (2 p числа).

На завершающем этапе каждый процессор одновременно и незави симо ведет расчет по формуле (3.20), для чего понадобится 4 m арифметических операций. В итоге получим оценку Tp (17 m 18 p)ta 8 ptcomm и T1 8 p / Sp.

T p 1 18 p / 17 n 8 appa * p 2 / 17 n Здесь ta max(tadd, tmult, tdiv ) – характерное время выполнения од ной арифметической операции с плавающей точкой;

tcomm – время на передачу одного числа между процессорными элементами, appa tcomm / ta.

Рис. 3.14. Оценка ускорения метода параметрической прогонки На рис. 3.14 приведены результаты расчетов ускорения по приве денной выше формуле. Из рисунка видно, что при небольших раз мерах задачи и низких показателях пропускной способности комму никационной среды, связывающей узлы многопроцессорной вычис лительной системы с распределенной памятью, использование рас смотренного метода параллельного решения системы (3.13) не явля ется эффективным.

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

Метод циклической редукции При решении системы линейных алгебраических уравнений с трехдиагональной матрицей методом циклической редукции требу ется, чтобы n 2q, q. Причём если первое и последнее уравнения в (3.11) имеют вид b0 x0 c0 x1 f0 ;

an xn-1 bn xn f n и c0 0 и an 0, то для удобства реализации метода циклической редукции рекомендуется расширить систему путем добавления до полнительных уравнений x0 0, xn 0 со сдвигом индексации со ответствующих неизвестных.

Итак, пусть система трёхдиагональных уравнений (3.13) имеет вид x0 f0, q ai xi 1 bi xi ci xi 1 fi, i 1,..., n 1;

n 2 ;

(3.25) x f.

n n Идея метода циклической редукции для решения этой задачи за ключается в следующем. Циклически, на каждом шаге k 1,..., q производится с помощью линейных преобразований исключение не известных из системы (3.25) таким образом, что остаются неизвест ные с индексами, кратными 2k. Например, на первом шаге ( k 1 ) из системы исключаются неизвестные с нечетными индексами следую щим образом ( i 2, 4,..., n 2 ):

ai 1 xi 2 bi 1 xi 1 ci 1 xi fi 1, (3.26) ai xi 1 bi xi ci xi 1 f i, ai 1 xi bi 1 xi 1 ci 1 xi 2 fi 1.

ai c, а третье – на i и при Умножим первое уравнение (3.26) на bi 1 bi бавим их ко второму уравнению. В результате получим следующую систему n / 2 1 уравнений, содержащую неизвестные с четными индексами:

x0 f0, (1) (1) (1) (1) (3.27) ai xi 2 bi xi ci xi 2 fi, i 2, 4,..., n 2, x f.

n n ai 1ai (1) ca ac cc ai(1), bi bi i 1 i i 1 i, ci(1) i 1 i, Здесь bi 1 bi 1 bi 1 bi fa fc fi (1) fi i 1 i i 1 i ;

i 2, 4,..., n 2.

bi 1 bi Продолжая редукцию неизвестных, на последнем (q 1) -м шаге по лучим систему, содержащую только три уравнения:

x0 f0, ( q 1) ( q 1) ( q 1) ( q 1) (3.28) an /2 x0 bn /2 xn / 2 cn /2 xn f n /2, x f.

n n Найденное из (3.28) значение f n(/21) anq 1) x0 cnq2 1) xn q ( ( /2 / xn / bnq 1) ( / может затем использоваться для определения значений неизвестных xn /4, x3 n /4 в редуцированной на ( q 2) -м шаге системе x0 f0, ( q 2) ( q 2) ( q 2) ( q 2) an / 4 x0 bn /4 xn /4 cn /4 xn /2 f n /4, ( q 2) ( q 2) ( q 2) ( q 2) (3.29) an / 2 xn /4 bn /2 xn /2 cn /2 x3 n /4 f n /2, ( q 2) ( q 2) ( q 2) ( q 2) a3 n /4 xn / 2 b3n /4 x3n /4 c3 n /4 xn f3n /4, xn f n.

Пользуясь таким приемом, можно последовательно найти остальные неизвестные с индексами, кратными n / 8, n / 16,..., 2,1.

Следовательно, процедура циклической редукции включает вы числения для шагов k 1,..., q 1 по рекуррентным формулам новых коэффициентов и правых частей:

ai( k ) Pai( 2k1)1 ;

bi( k ) bi( k 1) Pci( 2k1)1 Qi ai( 2k1)1 ;

ci( k ) Qi ci( 2k1)1 ;

k k k k i i ai( k 1) c ( k 1) fi ( k ) f i ( k 1) Pi fi (k21) Qi f i k21) ;

( ;

Qi i( k 1) ;

Pi (3.30) k 1 k bi(k2k1)1 bi 2k для i 2k,2 2k,..., n 2k при ai(0) ai, bi(0) bi, ci(0) ci, fi (0) fi с последующим замещением решения для k q, q 1,..., 2,1 из fi ( k 1) ai( k 1) xi 2k 1 ci( k 1) xi 2k xi, (3.31) bi( k 1) где i 2k 1,2 2k 1,..., n 2k 1.

На рис. 3.15 представлена диаграмма маршрутизации прямого и обратного хода циклической редукции для n = 8. Оценим временные затраты данного алгоритма.

Для расчета всех коэффициентов каждого уравнения редуциро ванной системы потребуется 12 арифметических операций, всего таких расчетов необходимо сделать (n / 2 1) ( n / 4 1)... 1 n q 1. При вычислении значения од ного неизвестного используется 5 арифметических операций, а всего их будет 5n. В итоге метод циклической редукции, также как и метод прогонки, является экономичным методом*, которому требу ется около 17n арифметических операций.


Рис. 3.15. Диаграмма маршрутизации прямого и обратного хода метода циклической редукции Этот алгоритм обладает хорошим параллелизмом (рис. 3.16). Во первых, операции исключения (3.30) независимы и могут выполнять ся параллельно. Во-вторых, расчет значений неизвестных по форму лам (3.31) проводится одновременно. Однако если используется мно гопроцессорная система с распределенной памятью и обычная схема декомпозиции данных по процессорам, то на этапе редукции понадо бится передача полученных коэффициентов в локальную память других процессорных элементов. Кроме того, на завершающих шагах прямого хода будет удваиваться число простаивающих процессоров.

Похожая картина наблюдается и при расчете значений неизвестных.

Сначала вычисляется значение xn /2 на одном процессоре, затем пе ресылается остальным, далее одновременно на разных процессорных элементах рассчитываются xn /4, x3 n /4 и т.д.

* Численный метод называется экономичным, если он корректен и число арифметических операций в нем пропорционально размеру задачи.

Рис. 3.16. Параллельная реализация метода циклической редукции Оценим трудоёмкость полученной параллельной реализации ал горитма метода циклической редукции для случая p 2r n. Пусть при инициализации в локальной памяти каждого процессорного эле мента размещены по m n / p уравнений системы (3.25) (за исклю чением первого или последнего, где на одно уравнение больше). То гда на этапе редукции (прямом ходе) временные затраты на вычисле ния коэффициентов оцениваются как 12 m ta, к этому нужно доба вить время на выполняющиеся одновременно на каждом шаге редук ции пересылки набора коэффициентов Fi 2k 1 ai 2k 1, bi 2k 1, ci 2k 1, fi 2k 1 между процессорными элементами, номера которых отличаются на единицу. С учетом оценки времени выполнения параллельного алгоритма на этапе вычисления неиз вестных получим Tp 17m ta 5q tcomm и ускорение t T1 p, appa comm.

Sp T p 1 5 appa * p log 2 ( n) ta 17 n Рис. 3.17. Оценка ускорения метода циклической редукции На рис. 3.17 представлены графики оценки ускорения параллель ного алгоритма циклической редукции. Сравнивая рис. 3.14 и 3.17, можно сделать вывод, что параллельный алгоритм метода цик лической редукции при небольших размерах задачи ( n ~ 102 ) и для многопроцессорных систем с медленной скоростью передачи данных между вычислительными узлами ( appa tcomm / ta ~ 103 ) так же, как и алгоритм параметрической прогонки, не является эффективным.

При увеличении размера задачи до n ~ 103 рассмотренный парал лельный вариант метода циклической редукции предпочтительнее, поскольку позволяет получить более высокие оценки величины ус корения параллельной программы.

Рассмотрим вариант метода циклической редукции, поддержи вающий параллелизм операций на максимально высоком уровне, что обеспечивается за счет другого способа организации вычислений на фазе редукции. На первом шаге из уравнений системы поочередно исключаются неизвестные с четными и нечетными номерами. В ре зультате получаются две независимые подсистемы. На втором шаге редуцируется похожим образом каждая полученная подсистема, что приводит в результате к четырем независимым подсистемам, содер жащим по (n 1) / 4 неизвестных. На последнем шаге этапа система уравнений (3.25) преобразуется к системе с диагональной матрицей (или к n 1 независимым уравнениям для отдельных неизвестных), решение которой можно получить, проводя одновременно вычисле ния по явным формулам.

Рис. 3.18. Диаграмма маршрутизации метода циклической редукции без обратного хода [1] На рис. 3.18 приведена диаграмма маршрутизации построенного алгоритма, который часто называют методом циклической редукции без обратного хода. Из диаграммы видно, что рассмотренный метод обладает максимальным параллелизмом, что стало возможным за счет увеличения общего числа арифметических операций. Следует также отметить, что одновременно существенно возросли и комму никационные затраты на этапе редукции, связанные с более активной пересылкой коэффициентов и правых частей линейных уравнений.

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

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

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

4. ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ ЛИНЕЙНЫХ СИСТЕМ Итерационные методы для решения линейных систем Ax b с невырожденной квадратной матрицей A nn используя заданное начальное приближение x0 n выполняют его последовательное улучшение до тех пор, пока приближенное решение не будет найде но с требуемой точностью. Итерации заканчиваются, когда норма невязки n b Axk max bi aij x (j k ) 1i n j или другая мера ошибки (например, относительная или абсолютная погрешность определения компонент решения ) не станет малой.

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

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

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

- линейную комбинацию векторов (saxpy);

- скалярное произведение векторов (dot);

- умножение матрицы на вектор (gaxpy);

- решение систем с треугольными матрицами.

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

мы cond ( A) A A параллельных процессов, которая сопровождается передачей данных между процессорными элементами для обеспечения локальных вы числений.

4.1. Метод Якоби Задавая начальное приближение x0, в методе Якоби новое ( k 1 )-е приближение к точному решению для каждой компоненты рассчитывается по следующей формуле ( aii 0, i 1,2,..., n ):

n bi aij x (jk ) j j i xi( k 1), i 1,..., n;

k 0,1,.... (4.1) aii Здесь i – номер компоненты вектора приближенного решения xk ;

k – номер итерации.

Пусть D – диагональная матрица, образованная диагональными элементами A, а L и U – нижняя и верхняя треугольные матрицы вида 0 0 0 a12 a1n a21 0 0, U 0 0 a2 n, D L U A, L an1 an 2 0 0 0 тогда (4.1) можно записать в виде xk 1 D 1 L U xk D 1b. (4.2) Для корректности метода Якоби требуются ненулевые диагональ ные элементы матрицы A, что можно обеспечить перестановкой строк или столбцов такой матрицы, если это необходимо. При реали зации метода Якоби на компьютере в памяти должны храниться все компоненты векторов xk и xk 1, поскольку ни один из компонентов xk не может быть переписан до тех пор, пока новое приближение xk 1 не будет получено. Кроме того, следует обратить внимание, что по формулам (4.1) или (4.2) все компоненты следующего приближе ния xk 1 могут быть рассчитаны одновременно и независимо, что открывает большие возможности при построении параллельной вер сии такого алгоритма.

Конечно, метод Якоби сходится не всегда, но если матрица A имеет строгое диагональное преобладание n n aii aij, i 1,..., n или a jj aij, j 1,..., n, j 1 i j i i j то этого достаточно для сходимости метода, хотя сходимость может быть очень медленной. Кроме того, для проверки сходимости метода Якоби используются необходимые и достаточные условия: корни характеристического уравнения a11 a12 a1n a21 a22 a2 n det an 2 ann an должны по модулю быть меньше единицы.

В качестве критерия завершения итерационного процесса рас сматривают следующее условие:

xk 1 xk max xi( k 1) xi( k ) и Axk 1 b, (4.3) i где xk – приближенное значение решения на k -й итерации числен ного метода.

Пользуясь формулой (4.1), оценим временные затраты на выпол нение одной итерации рассмотренного алгоритма метода Якоби (4.1):

T1 n ( n 1) tadd tmult t add tdiv 2 n 2 ta ;

(n 1), где ta – обобщенное время одной арифметической операции (+, –, *, /).

При построении параллельной версии итерационного метода Яко би будем пользоваться подходом, опирающимся на декомпозицию данных, т.е. для p-процессорной параллельной системы будем на значать на расчет (k+1)-го приближения по формуле (4.1) по n/p компонентов вектора xk 1, для чего понадобится n/p строк матрицы A и n/p компонентов вектора b. К такому же способу организации па раллельных вычислений можно прийти, если пользоваться стандарт ной технологией разработки параллельных программ, подробно опи санной во введении.

Действительно, на этапе декомпозиции в качестве фундаменталь ной мелкозернистой подзадачи выберем задачу расчета одного ком понента вектора xk 1 по известным значениям вектора предыдущего приближения xk, коэффициентов матрицы A и вектора правой части b (4.1).

При проектировании коммуникаций между подзадачами стано вится очевидным, что они не требуются, поскольку расчеты каждого компонента xk 1 на одной итерации ведутся независимо и одновре менно. Однако следует отметить, что при выполнении любого итера ционного метода необходимо контролировать выполнение условия (4.3), а для этого (для вычисления нормы ошибки и нормы невязки) нужно располагать всеми компонентами полученного нового при ближения xk 1, а также значениями xk, A, b. Таким образом, потре буется сборка xk 1 либо на одном или на всех процессорных элемен тах с последующим расчетом условия (4.3) и принятия решения о продолжении итерационного процесса.

Первый способ (сборка xk 1 на одном процессорном элементе) использует меньшее количество межпроцессорных коммуникаций, однако, когда будут рассчитываться значения норм в (4.3), остальные процессоры будут простаивать, что снизит эффективность парал лельной программы. Кроме того, после вычисления логического зна чения по неравенствам (4.3) процессорный элемент, который прово дил вычисления, должен передать сообщение остальным процессор ным элементам системы о принятом решении о продолжении или остановке итерационного процесса, на что также потребуются об менные операции.

При сборке xk 1 на каждом процессорном элементе потребуется несколько большее количество межпроцессорных коммуникаций, но многие из них при топологии многопроцессорной вычислительной системы с высокой степенью связности могут проходить одновре менно (здесь также требуется, чтобы на каждом ПЭ были xk, A, b ).

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

Кроме того, поскольку для продолжения вычислений на следующей итерации требуется иметь текущее приближение, то этот способ име ет преимущество.

Таким образом, хотя выбранный способ декомпозиции обеспечи вает получение нового приближения к точному решению системы Ax b без межпроцессорных коммуникаций, проверка условия за вершения итерационного процесса (4.3) является узким местом раз рабатываемого параллельного алгоритма/программы. Это связано с наличием глобальных коммуникаций в массиве фундаментальных подзадач, которые необходимы для расчета условия (4.3) и для об новления вектора последнего приближения xk в (4.1). Также заме тим, что требование хранения всей матрицы A на одном или на каж дом ПЭ ограничивает размер задачи n, который может быть значи тельным.

Рассмотрим схему укрупнения мелкозернистых фундаментальных подзадач, при которой каждый укрупненный блок содержит по n/p таких подзадач. Для такой схемы укрупнения расчет n/p компонентов xk 1 в каждом блоке не требует хранения всей матрицы A, а лишь n/p ее строк (см. рис. 4.1). Это же касается вектора b – для расчетов нужны лишь n/p его компонентов. Вектор xk должен быть представ лен всеми своими компонентами (рис. 4.1).

Параллельная версия итерационного метода Якоби в таком случае реализуется следующим образом:

- Этап инициализации: на каждом процессорном элементе ПЭ размещается A, b ( 1,..., p ) и x0, где A, b содержат стро ки матрицы A и компоненты вектора b с индексами n n i 1 ( 1),..., ;

1,..., p.

p p - Этап параллельных вычислений:

для k 0,1, 2,... выполнять на каждом ПЭ:

{1} расчет xk 1 по формуле (4.1);

{2} расчет xk 1 xk и b A xk ;

{3} пересылка полученных значений норм вместе с вектора ми xk 1 на остальные ПЭ;

{4} обновление вектора последнего приближения и проверка каждым ПЭ выполнения условия (4.3). В случае, если условие выполнено, то выход из цикла.

Ap, b, A1, b, A2, b, p 1 xk, xk 1 1 xk, xk 1 2 xk, xk 1 p … ПЭ 1 ПЭ 2 ПЭ p Рис. 4.1. Распределение исходных данных (k = 0) по процессорным элементам Оценим временные затраты разработанного параллельного алго ритма на одной итерации без учета этапа инициализации n / p 1 :

n2 n n n n ta 4 ta p 1 2 tcomm 2 ta p 1 tcomm Tp.

p p p p p {1} {2} {3} При оценке коммуникационных затрат предполагается, что каж дый ПЭ одновременно рассылает по n / p 2 числа остальным про цессорным элементам.

Тогда получим оценку ускорения и эффективности алгоритма:

2 n 2t a T1 p (4.4) Sp ;

( p 1) Tp n n 2 ta p 1 tcomm 2n p p Sp t, comm.

Ep p 1 ( p 1) ta 2n Из формул (4.4) следует, что можно достичь высокой эффектив ности распараллеливания, если n / p 1, tcomm / ta – показа тель, характеризующий соотношение быстродействия передачи дан ных по межпроцессорной сети к быстродействию процессоров. Сле дует отметить еще одно свойство итерационного метода Якоби: по следовательная и параллельная версии алгоритма выполняют одина ковый объем вычислительной работы, что характеризуется одинако вым количеством итераций для получения приближенного решения с заданной точностью.

Рис. 4.2. Ускорение параллельного метода Якоби при решении СЛАУ с полнозаполненной матрицей (n = 5000). Линия – оценка по формуле (4.4) при = 140. Значки – время расчета по параллельной программе на кластере СКИФ Cyberia На рис. 4.2 представлены результаты зависимости оценки ускоре ния описанного выше параллельного алгоритма и построенной по нему параллельной программы для n = 5000 и = 510-6 от числа ис пользуемых процессоров p. Численные расчеты проведены на супер компьютере СКИФ Cyberia при решении системы линейных уравне ний с матрицей, все элементы которой равны единице за исключени ем диагональных, имеющих значение 1.1*(n – 1).

Для полученных по формуле (4.4) теоретических оценок ускоре ния параллельного алгоритма = 140. Из рисунка видно, что рост ускорения параллельной программы с увеличением числа исполь зуемых процессорных элементов p уменьшается, что связано (см.

левую часть (4.4)) с сокращением объема вычислений при неизмен ном объеме пересылаемых между процессорными элементами дан ных.

4.2. Метод Зейделя и метод верхней релаксации Более быструю сходимость при последовательных вычислениях имеет метод Зейделя. Оценки показывают, что скорость сходимости метода Зейделя в случае диагонального преобладания матрицы A в среднем в 2–3 раза выше, чем метода Якоби.

Покомпонентная форма записи метода Зейделя имеет вид ( aii 0, i 1,2,..., n ):

i 1 n bi aij x (jk 1) (k ) ax ij j j 1 j i xi( k 1), i 1, 2,..., n. (4.5) aii Однако это не единственная форма записи этого метода, основная идея которого заключается в использовании только что вычисленных на ( k 1 )-й итерации компонент вектора xk 1 для расчета остальных.

Например, можно рассмотреть и такую форму:

i 1 n bi aij x (jk ) ( k 1) a x ij j j 1 j i xi( k 1), i n, n 1,...,1. (4.6) aii Также можно рассматривать и другой порядок изменения значе ния индекса i, главное, чтобы только что вычисленные xi( k 1) сразу же использовались в вычислительном процессе.

Матрично-векторная форма записи (4.5) и (4.6) имеет вид:

xk 1 D L b Uxk, (4.7) xk 1 D U b Lxk, (4.8) где D, L, U – диагональная, нижнетреугольная и верхнетреугольная матрицы, составленные из элементов матрицы A, причем D L U A.

Метод Зейделя сходится при любом начальном приближении и любой правой части b, если матрица A имеет строгое диагональное преобладание. Необходимые и достаточные условия сходимости ме тода Зейделя заключаются в требовании, чтобы все корни характери стических уравнений a11 a12 a1n a21 a22 a2 n det 0 для метода (4.7) an1 an 2 ann или a11 a12 a1n a21 a22 a2 n det 0 для метода (4.8) ann an1 an были бы по модулю меньше единицы.

Анализируя формулы (4.5) и (4.6), можно заметить, что в отличие от метода Якоби для хранения xk 1 и xk требуется один вектор, и расчет по формулам метода Зейделя подразумевает лишь последова тельное выполнение входящих в формулы арифметических опера ций. Дело в том, что формулы (4.5) и (4.6) являются рекуррентными, т.е. для расчета xi( k 1) в (4.5) нужно предварительно на ( k 1 )-й ите рации рассчитать значения x1( k 1), x2k 1),..., xi(11).

( k Так, если использовать представленную на рис. 4.1 схему распре деления исходных данных по процессорным элементам, то каждое только что рассчитанное значение xi( k 1) в соответствии с (4.5) не медленно должно быть отправлено на остальные процессорные эле менты, которые вынуждены простаивать в ожидании этого значения, и, следовательно, расчет каждого xi( k 1) ведется одним ПЭ при про стаивающих остальных (последовательное выполнение алгоритма метода Зейделя (4.5)).

Для преодоления этой проблемы в практике параллельных вычис лений часто применяют так называемую асинхронную, или хаотиче скую релаксацию, когда при получении компонент вектора прибли женного решения на следующей итерации используют последние значения xk и xk 1, имеющиеся на каждом процессорном элементе (рис. 4.1), а передача вычисленных значений другим процессорным элементам осуществляется после завершения итерации. Этот подход может быть эффективным, однако порядок выполнения вычисли тельной работы в нем отличается от организации расчетов в ориги нальном методе Зейделя, что существенно усложняет анализ сходи мости полученного параллельного метода решения систем линейных алгебраических уравнений.



Pages:     | 1 || 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.