Управление потоками данных в параллельных алгоритмах вычислительной линейной алгебры

Недостатком программ, созданных с помощью данных языков программирования, является их недостаточная переносимость. Кроме того, компиляторы этих языков не всегда входят в стандартный комплект программного обеспечения, поставляемый с компьютером. 2. межпроцессных коммуникаций ( interprocess communication ) такие как разделяемая память ( shared memory ), семафоры ( semaphores ), очереди сообщений ( message queues ), сигналы ( signals ) удобно использовать при программировании в модели разделяемой памяти. При использовнии низкоуровневых средств межроцессных коммуникаций программисту кроме кодирования непосредственно вычислительного алгоритма, требуется самому создавать процедуры синхронизации процессов, что может быть источником дополнительных трудностей. 3. Parsytec , поставляеются со специализированной операционной системой Parix , обеспечивающей реализацию параллельных алгоритмов на программном уровне. ОС Parix предоставляет разработчику библиотеку системных функций, обеспечивающих синхронную и асинхронную передачу данных, получение информации о конфигурации вычислительного узла, на котором запущена программа и т.д.

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

Обычно в вычислительных системах, управляемых потоком данных, команды машинного уровня управляются доступностью данных, проходящих по дугам графа потока данных (ГПД). В данной системе используется принцип управления укрупненным потоком данных ( Large-Grain Data Flow ), в котором единица планирования вычислений — процесс — крупнее (возможно, намного крупнее), чем одна машинная команда. Для задания ГПД был разработан специальный язык DGL ( Dataflow Graph Language ). ГПД, записанный на этом языке, легко настраивается под конкретную многопроцессорную систему — число вершин и дуг графа может определяться через внешние параметры, которые вводятся пользователем, считываются из файла или задавются каким-либо другим способом. В качестве параметров можно использовать число процессоров в системе, степень загруженности системы и т.д. В состав системы FLOWer входят: · · DGL ; · Структура параллельной программы, написанной в системе FLOWer , мало отличается от структуры последовательной программы. Так, процессы записываются ввиде отдельных функций, которые считавают значения параметров из входных дуг ГПД и передают результат в выходные.

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

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

Основные понятия параллелелизма Определение.Степенью параллелизма численного алгоритма называется число его операций, которые можно выполнять параллельно.

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

Определение.Ускорением параллельного алгоритма называется отношение где T 1 - время выполнения алгоритма на одном процессоре, Т p - время выполнения алгоритма в системе из p процессоров.

Параллельный алгоритм называется масштабируемым, если при увеличении числа процессоров производительность также увеличивается. В идеальном случае Определение.Ускорением параллельного алгоритма по сравнению с наилучшим последовательным алгоритмом называется отношение где T’ 1 - время выполнения алгоритма на одном процессоре, Т p - время выполнения алгоритма в системе из p процессоров.

Определение.Эффективностью параллельного алгоритма называется величина Очевидно, что S p SYMBOL 163 f 'Symbol' s 12 Ј p , S’ p SYMBOL 163 f 'Symbol' s 12 Ј S p , E p SYMBOL 163 f 'Symbol' s 12 Ј . Если алгоритм достигает максимального ускорения ( S p = p ), то E p = 1. Одной из целей при конструировании параллельных алгоритмов является достижение по возможности большего ускорения. Но максимальное ускорение можно получить только для тривиальных задач.

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

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

Балансировка нагрузки может осуществляться как статически, так и динамически.

Рассмотрим формальную модель системы, в которой где SYMBOL 97 f 'Symbol' s 12 a 1 - доля операций, выполняемых только одним процессором, SYMBOL 97 f 'Symbol' s 12 a 2 - доля операций, выполняемых со средней степенью параллелизма k p , SYMBOL 97 f 'Symbol' s 12 a 3 - доля операций, производимых со степенью параллелизма p , t d - общее время, требуемое для подготовки данных.

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

Случай 1. SYMBOL 97 f 'Symbol' s 12 a 1 = 0, SYMBOL 97 f 'Symbol' s 12 a 2 = 0, SYMBOL 97 f 'Symbol' s 12 a 3 = 1, t d = 0. Здесь ускорение максимально.

Случай 2. SYMBOL 97 f 'Symbol' s 12 a 1 = 0, SYMBOL 97 f 'Symbol' s 12 a 2 = 1, SYMBOL 97 f 'Symbol' s 12 a 3 = 0, t d = 0. Теперь S p = k p . Случай 3. SYMBOL 97 f 'Symbol' s 12 a 1 = SYMBOL 97 f 'Symbol' s 12 a SYMBOL 97 f 'Symbol' s 12 a 2 = 0, SYMBOL 97 f 'Symbol' s 12 a 3 = 1 - SYMBOL 97 f 'Symbol' s 12 a t d = 0. В этом случае Данная формула выражает закон Амдаля-Уэра . Следует отметить очень быстрое убывание S p для малых значений SYMBOL 97 f 'Symbol' s 12 a Используя последнюю формулу, можно примерно оценить долю параллельных вычислений Случай 4. t d велико. Этот случай показывает, что при достаточно большом t d можно получить S p При вычислении ускорения по описанным формулам возникает ряд трудностей. В частности, задача, решаемая с большим числом процессоров может не помещаться в оперативной памяти единственного процессора; исследование ускорения должно проводится для больших задач, которые слишком велики для того, чтобы их можно было решить за разумное время с помощью единственного процессора.

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

Рассмотрим, например, задачу матрично-векторного умножения Ax . Если An*n-матрица , то потребуется 2n 2 - n операций.

Предположим, что при исходном порядке n вычисления на единственном процессоре идут со скоростью 1 MFLOP. Затем мы удваиваем n , и число операций при большом n увеличивается примерно в 4 раза. Если задача решается на 4 процессорах со скоростью 4 MFLOP, то мы достигли максимального ускорения.

Вообще, пусть размер задачи определяется параметров n , а число операций есть f ( n ) при растущем n и p = SYMBOL 97 f 'Symbol' s 12 a f ( n ). Если n 1 - размер задачи, посильный для одного процессора, то соответствующий размер задачи n p для системы из p процессоров определяется равенством f ( n p ) = pf (n 1 ) (чтобы число операций, выполняемых p процессорами, в p раз превосходило число операций одного процессора). Определение.Параллельный алгоритм обладает свойством локальности, если он использует локальные данные гораздо чаще, чем удаленные.

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

Важность этого свойства определяется отношением стоимостей удаленного и локального обращения к памяти. Это отношение может варьироваться от 10:1 до 1000:1 и больше, что зависит от используемой аппаратуры. Метод оценки производительности Задачей программиста является проектирование и реализация программ, которые удовлетворяют требованиям пользователя в смысле производительности и корректной работы.

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

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

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

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

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

Прежде всего данные загружаются в буферную память или собираются в последовательных ячейках памяти процессора Р 1 . Затем выполняется команда переслать, результатом которой является перенос данных в буфер или память процессора Р 2 . Далее процессор Р 2 выполняет команду принять, и данные направляются на места своего конечного назначения в памяти Р 2 . Чтобы освободить основные процессоры от значительной части этой работы, можно использовать сопроцессоры. В этом упрощенном описании опущены многие детали (в большей степени зависящие от конкретной системы), но главные пункты сохранены: чтобы переслат данные, их вначале нужно выбрать из памяти передающего процессора, должна быть предоставлена информация, куда следует их послать, должен произойти физический перенос данных от одного процессора к другому и, наконец, данные нужно разместить в надлежащих ячейках памяти принимающего процессора. Для многих систем время такого обмена можно выразить приближенной формулой где T start — время запуска, а T send — время, необходимое для пересылки одного числа.

Подобным образом оценивается время перемножения и сложения n пар чисел. Будем считать, что n пар чисел перемножаются за время nT mul , а складываются за время nT add . При этом не будем различать операции сложения и вычитания, и операции умножения и деления.

Экспериментально было установлено, что время T mul на порядок превосходит время T add (конкретное значение зависит от типа процессора и отношение T mul / T add варьируется примерно от 10 до 100). Поэтому часто при оценке времени выполнения алгоритма операции сложения не учитываются. Ч А С Т Ь 1 СИСТЕМА FLOWer В данной части описывается система программирования FLOWer — набор инструментов, в значительной степени облегчающих создание и отладку параллельных программ. Эта система написана на языке Delphi и предназначена для работы в локальных сетях ЭВМ под управлением ОС Windows . Г Л А В А 1 КРАТКИЙ ОБЗОР Программный комплекс FLOWer спроектирован для работы в массивно-параллельной системе ( MPP ), которая состоит из однородных вычислительных узлов, включающих в себя: · · · Узлы связаны через некоторую коммуникационную среду (высокоскоростная сеть, коммутатор и т.п.). На каждом зле работает полноценная ОС (в данном случае Windows ). В частности, такой системой является локальная сеть, в которой любые два компьютера могут обмениваться данными друг с другом. На любом компьютере могут одновременно работать несколько процессов.

Данная работа была выполнена в локальной сети с общей шиной. В такой сети все процессоры соединены посредством шины.

Достоинством такой системы является очень малое число линий связи, однако возможны задержки (конфликт на шине) при использовании шины несколькими процессорами. Это может стать серьезной проблемой при увеличении числа процессоров. В состав системы FLOWer входят: · · DGL ; · Запуск параллельной программы осуществляется под управлением диспетчера, который распределяет процессы по компьютерам и устанавливает связи между процессами. Для нормальной работы диспетчера на всех компьютерах должна быть запущена специальная программа - монитор.

Создание параллельной программы в системе FLOWer состоит из следующих шагов: · · DGL · DGL · · DLL · DGL Г Л А В А 2 МОДЕЛЬ ВЫЧИСЛЕНИЙ Система FLOWer базируется на модели управления потоком данных. В данной модели вычислений программа представляется в виде ГПД. Вершины ГПД соответствуют отдельным процессам — последовательным участкам программы, а дуги задают отношения между ними. Точка вершины, в которую входит дуга, называется входным портом (портом импорта или входом), а точка, из которой она выходит, - выходным (портом экспорта или выходом). Некоторые входные порты могут быть помечены как обязательные. По дугам передаются данные из одного процесса в другой.

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

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

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

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

Каждому процессу сопоставлена пара чисел, однозначно определяющая этот процесс, — номер класса процесса и номер экземпляра процесса (о необходимости введения двумерной нумерации процессов будет рассказано несколько позже). Различие между классом и экземпляром процесса такое же, как и в объектно-ориентированных языках программирования. Канал определяет однонаправленный поток данных между процессами, он всегда направлен от входа канала к выходу. Между каналами и выходами имеется взаимооднозначное соответствие; в то же время вход может быть соединен с несколькими выходами.

Каждому каналу приписывается вес — величина, характеризующая примерный объем данных, пересылаемых по каналу.

Процесс — структура данных, состоящая из набора входов и выходных портов (не путать с выходами). Входы и выходные порты занумерованы. Вход может быть помечен как обязательный.

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

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

Введем систему обозначений. · · Pi — процессы с номером ( i , x ), где x Pi |. · Pij — процесс с номером ( i , j ). · In Pij — множество входов процесса ( i , j ). · In Pij | — число входов процесса ( i , j ). · In k Pij — k-й вход процесса ( i , j ). · Out Pij — множество выходов процесса ( i , j ). · Out k Pij — выходы процесса ( i , j ) с номерами ( k , x ), где x Out k Pij | · Out kl Pij — выход с номером ( k , l ). · Out kl Pij , тогда выход X соединен с входом h = X.in процесса с номером ( s , w ), где s = X.proc , w = X.copy . Т.е. существует канал c = (X, Y), где Y = In h Psw . 2.2. Шаблон ГПД Шаблоны ГПД имеют сходные черты с шаблонами в языке С++. Каждый шаблон ГПД обладает набором параметров, после определения которых из шаблона получается конкретный экземпляр ГПД. От параметров может зависеть число процессов ГПД, число выходов процесса.

Параметры шаблона вводятся на этапе запуска программы. Это позволяет настраивать программу под конкретные требования пользователя. Для записи шаблона ГПД был разработан специальный язык DGL ( Dataflow Graph Language ). Более подробно о DGL можно прочитать в главе «Язык DGL ». Шаблон ГПД — это ориентированный граф, вершинами которого являются классы процессов, а дугами — каналы, и набор параметров.

Каждому классу процесса сопоставлен номер.

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

Введем систему обозначений. · · Ai — i-й параметр (целое число). · TPi — i-й процесс. · In TPi — множество входов процесса i . · In k TPi — k-й вход процесса i . · Out TPi — множество выходов процесса i . · Out k TPi — k-й выход процесса i . Пусть X = TPi . Тогда · X.ncopies = ncopies (A) — целочисленная функция от параметров, которая задает число копий процесса X. Пусть X = Out k TPi . Тогда · X.ncopies = ncopies ( p , A) — целочисленная функция, которая задает число копий выхода X. Здесь p — целое число, А — параметры. · X.proc — номер процесса, с входом которого соединен данный выход. · X.copy = copy ( c , p , A) — целочисленная функция. Здесь p и c — целые числа, А — параметры. · X.in — номер входа, с которым соединен данный выход. 2.3. Связь ГПД и шаблона ГПД Пусть TG (TP, TC, TA) — шаблон ГПД, A — конкретные значения параметров.

Введем операцию подстановки Subst параметров А в шаблон TG. G = Subst A TG — ГПД такой, что · TPi сопоставлено множество процессов Pi . · Pi | = ( TPi ). ncopies (A). · In Pij = In TPi , где j · Out k TPi сопоставлено множество выходов Out k Pij , где j Pi |. · | Out k Pij | = (Out k TPi ). ncopies (j, A), где j · (Out kl Pij ).proc = (Out k TPi ).proc, где l k Pij |, j · (Out kl Pij ).copy = (Out k TPi ).copy (l, j, A), где l k Pij |, j · (Out kl Pij ).in = (Out k TPi ).in, где l k Pij |, j ЗАМЕЧАНИЕ: при некоторых значениях вектора параметров операция подстановки может быть не определена. Г Л А В А 3 ЯЗЫК DGL Язык DGL был специально разработан для описания ГПД. Для ввода програмы на языке DGL можно использовать как простой текстовый редактор, так и специальную программу Rose , которая рисует ГПД на экране компьютера.

Рассмотрим как зписываются графы на языке DGL на основе элементарных примеров. 1. Out вершины V 1 связан со входом In вершины V 2 . process V1 { export: Out V2:In; } process V2 { import: In; } 2. Имеется N вершин V 0 , …, V N-1 . Вершина W содержит N выходов Out 0 , …, Out N -1 . Выход Out k соединен со входом In вершины V k , k = 0, …,( N -1). process V[N] { import: In; } process W { export: Out [N] V[c]:In; } 3. Имеется N вершин V 0 , …, V N -1 . Выход Out вершины k соединен со входом In вершины W , k = 0, …,( N -1). process W { import: In; } process V[N] { export: Out W:In ; } 4. Имеется N вершин V 0 , …, V N -1 и N вершин W 0 , …, W N -1 . Выход Out вершины V k связан со входом In вершины W k , k = 0, …, ( N -1). process V[N] { export: Out W[p]:In; } process W[N] { import: In; } 5. Имеется N вершин V 0 , …, V N -1 и M вершин W 0 , …, W M -1 . Каждая вершина типа V содержит M выходов Out 0 , …, Out M -1 , причем выход Out k связан со входом In вершины W k , k = 0, …, ( M -1). process W[M] { import: In; } process V[N] { export: Out [M] W[c]:In; } 6. Имеется N вершин V 0 , …, V N -1 . Выход Out вершины V k связан со входом In вершины V ( k +1) mod N , k = 0, …, ( N -1). process V[N] { export: Out V[(p+1) mod N]:In; import: In; } Значение N может быть объявлено как внешний параметр. В этом случае N определяется во время запуска программы. В дальнейшем планируется ввести в язык многомерные массивы процессов, а так же добавить условные операторы if … then … else и case … of . Это позволит описывать более сложные ГПД, в частности, многомерные решетки.

Синтаксис DGL ::= «DATAFLOW» «PROGRAM» name «;» {«EXTERN» name [«=» integer] «;»} {«CONST» name «=» expression «;»} {process} process ::= «PROCESS» name [«[» number_of_instances «]»] {directive «;»} «{» { process_attr «;»} «EXPORT:» {export} «IMPORT:» {import} «}» directive ::= «LOCAL» | «START» | «TERMINATION» process_attr ::= «weight» «=» integer import ::= name [«{» import_attr {«;» import_attr } «}»] «;» export ::= name [«[» number_of_instances «]»] « --> » process_name [«[» process_index «]»] «:» import_name [«{» export_attr {«;» export_attr } «}»] «;» number_of_instances ::= expression process_index ::= expression import_attr ::= «ARGUMENT» export_attr ::= «DATASIZE» «=» integer Г Л А В А 4 ПРИМЕР ПАРАЛЛЕЛЬНОЙ ПРОГРАММЫ В качестве примера рассмотрим задачу приближенного вычисления числа Пи с использованием правила прямоугольников для вычисления определенного интеграла где Согласно правилу прямоугольников, где Следует отметить, что это «процессорная» программа. Она не затрагивает многие проблемы параллельного программирования, например критическое влияние процессов ввода-вывода. Тем не менее эта задача поможет ознакомиться с основными принципами построения программ, работающих в соответствии с методом управления потоком данных.

Существует множество подходов к решению контрольной задачи.

Решение, приведенное ниже, иллюстрирует все основные шаги разработки программы. 4.1. Конструирование ГПД программы В пределе разрабатываемая программа может быть создана в виде одного процесса, но при этом теряется параллелелизм . Можно создать множество мелких процессов, таких как один оператор или даже одна арифметическая операция, что приведет к резкому увеличению расходов, связанных с запуском каждого процесса и обменом данных между ними.

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

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

Каждой вершине ГПД соответствует процесс DGL с тем же именем.

Теперь нарисуем ГПД с помощью программы Rose и запишем его на диск в формате DGL . DATAFLOW PROGRAM PI; EXTERN NProcs = 5; CONST nw = MAX (1, NProcs - 2); PROCESS Summer [1]; LOCAL; { IMPORT: NumIter {}; PartSum {}; } PROCESS Worker [ nw ]; { EXPORT: PartSum [1] --> Summator [0] : PartSum {DATASIZE = 0}; Demand [1] --> Manager [0] : DemandList {DATASIZE = 0}; IMPORT: Arg {}; } PROCESS Manager [1]; LOCAL; START; TERMINATION; { EXPORT: Done [1] --> Summator [0] : DoneSignal {DATASIZE = 0}; Worker [ nw ] --> Worker [c] : Arg {DATASIZE = 0}; IMPORT: DemandList {}; } 4.3. Компиляция процессов После того, как ГПД записан на языке DGL , его нужно обработать транслятором dglc . В результате работы транслятора на диске должны появиться следующие файлы: PI . dpr , Manager . pas , ManagerBody . pas , Worker . pas , WorkerBody . pas , Summator . pas , SummatorBody . pas . Далее, программист записывает код каждого процесса в файлах SummatorBody . pas , WorkerBody . pas , ManagerBody . pas , а затем компилирует проектный файл PI . dpr . Ч А С Т Ь 2 РЕАЛИЗАЦИЯ НЕКОТОРЫХ АЛГОРИТМОВ ВЛА В СИСТЕМЕ FLOWer В данной части рассматривается реализация в системе программирования FLOWer некоторых алгоритмов ВЛА, в число которых входят умножение матриц, LU -разложение, решение треугольных СЛУ, QR -разложение, итерационные методы решения СЛУ. Для каждого алгоритма приводятся теоретические оценки ускорения, описание параллельного алгоритма, экспериментальные данные. Для проведения экспериментов использовалась локальная сеть FastEthernet , состоящая из 4-х компьютеров Pentium II 333 MHz . Г Л А В А 1 УМНОЖЕНИЕ МАТРИЦ Рассмотрим задачу вычисления произведений A x и AB , где А и В – матрицы, а х — вектор.

Отдельно будут рассмотрены ленточные и блочно-диагональные матрицы. 1.1. Умножение матрицы на вектор Пусть А – матрица m n , а х — вектор длины n .Тогда произведение можно записать двумя способами или где а i — i -я строка матрицы А, а i — i - й столбец матрицы А, а ( x , y ) — скалярное произведение.

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

Рассмотрим сначала векторно-матричное произведение с помощью линейных комбинаций: p = n и x i и a i приписаны процессору i .Все произведения x i a i вычисляются с максимальным параллелизмом, а затем выполняются сложения по методу сдваивания.

Причем синхронизация в данной модели не требуется.

Алгоритм скалярных произведений в некоторых отношениях более привлекателен. Пусть p = m , x и a i приписаны процессору i . Каждый процессор выполняет свое скалярное произведение, и при этом паралелизм максимален. Выбор алгоритма вычислений зависит от целого ряда обстоятельств.

Матрично-векторное умножение неизбежно является частью более широкого процесса вычислений, и основную роль играет способ хранения матрицы А и вектора х. Еще одно существенное соображение — желаемое расположение результата по окончании умножения: в первом случае результат размещается в памяти одного процессора, тогда как в другом он размещен между процессорами. В реальных системах n и m значительно больше числа процессоров, и каждому процессору передаются несколько строк или столбцов. 1.2. Матричное умножение Аналогично рассматриваются алгоритмы умножения матриц А и В. Пусть матрицы разбиты на блоки Пусть число процессоров р равно числу st блоков матрицы С. Тогда все блоки можно вычислять одновременно.

Рассмотрим специальные случаи разбиения: · s = 1 — матрица А разбита на группы столбцов; · t = 1 — матрица В разбита на группы строк; · s = t = 1 — алгоритм блочного скалярного произведения; · r = 1 — алгоритм блочного внешнего произведения. Для блочного внешнего произведения распределение блоков по процессорам будет выглядеть следующим образом: Для простоты будем считать, что матрицы А и В имеют размерность n n , n = sm , где s — число блоков, m — число строк (столбцов) в блоке. Далее, пусть умножение двух чисел происходит за время T mul , а пересылка одного числа — за время T send . Процесс вычисления состоит из трех этапов: · · вычисление произведения блоков ; · пересылка результата . Пусть t 1 ( n ) — время умножения матриц на одном процессоре, t p ( n ) — время умножения матриц в системе из p = s 2 процессоров. В случае последовательного алгоритма произведение матриц требует n 3 операций умножения, пересылки данных не требуются. В случае параллельного алгоритма пересылка исходных данных выполняется за время 2 nmT send , умножение блоков — за время m 2 nT mul , пересылка результата — за время n 2 T send . Тогда Построим график S p ( n ) при p =4 и при различных значениях отношения T send / T mul . При достаточно больших значениях n Т.е. получаем линейный прирост производительности при n ® . Совсем другая картина получается, если строить строить график S p ( n ) как функции от р при некотором фиксированном значении n . В этом случае S p ( n ) сначала резко возрастает до некоторого экстремального значения S ’= S p ’ ( n ), а затем начинает плавно убывать. Т.е. мы достигли некоторого оптимального значения числа процессоров p ’ и дальнейшее увеличение числа процессоров будет приводить только к замедлению программы. Это объясняется тем, что при увеличении числа процессоров увеличивается и объем передаваемых данных.

Значение p ’ в данном конкретном случае несложно вычислить. Оно находится из уравнения p ³ 1 Решая уравнение, получим p = p ’ максимум эффективности E p скорее всего не будет достигнут. Результаты эксперимента

Очевидно, что А n можно вычислить следующим образом Тогда вполне естественно приписать каждый блок отдельному процессору. При этом, если размерности блоков сильно отличаются, возникает проблема балансировки нагрузки.Данную проблему можно разрешить при помощи динамического распределения блоков между процессами. Пусть все А ii — матрицы размерности m m , тогда в системе из p = s процессоров Положим p =4, m =10 и построим график S p ( n ) при различных значениях отношения T send / T mul . Видно, что рост производительности параллельной системы наблюдается как при больших значениях m , так и при больших значениях n . Вычислим оптимальное число процессоров p ’ и построим график S p ( n , m ) при n =5000, m =10. Значение p ’ определяется из уравнения , p ³ 1 Отсюда Результаты эксперимента

Заметим, что если полуширина ленты для А равна a , а для В равна b , то полуширина ленты для АВ в общем случае будет равна a + b . Пусть А — матрица размерности n n с полушириной ленты a , В — матрица размерности n n с полушириной ленты a . Причем А разбита на блоков Число ненулевых элементов матрицы B не превышает (2 a +1) n , в блоке A i — (2 a +1) m . В результирующей матрице будет не более (4 a +1) n ненулевых элементов. В системе из p процессоров построим вычисление по следующей схеме: первый процессор вычисляет А 1 В, второй — А 2 В и т.д.

Соответственно, матрицу А выгодно хранить по строкам, а матрицу В — по столбцам. В результате умножения получим матрицу, которая будет храниться в памяти по строкам. На вычисление одного ненулевого элемента матрицы АВ затрачивается время не более (2 a + 1) T mul . Тогда Видно, что ускорение S p зависит от полуширины ленты a и не зависит от n . Теперь вычислим оптимальное значение числа процессоров p ’ и построим график функции S p ( a ) при a =500. Значение p ’ определяется из уравнения Отсюда Г Л А В А 2 ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ ЛИНЕЙНЫХ СИСТЕМ Рассмотрим систему линейных уравнений A x = b с невырожденной матрицей А размера n n . 2.1. LU -разложение Постановка задачи Построим разложение мтрицы A = LU , где L — нижнетреугольная матрица с единицами на главной диагонали, а U — верхнетреугольная матрица. В качестве алгоритма LU -разложения выберем kij -форму с выбором главного элемента по строке (столбцу). Упрощенно алгоритм можно записать следующим образом для k = 1 до n -1 для i = k+1 до n l ik = a ik /a kk для j = k+1 до n a ij = a ij - l ik a kj Распараллеливание и оценки производительности Предположим вначале, что система состоит из p = n процессоров. Тогда один из возможных вариантов организации LU -разложения выглядит так. Пусть i - ая строка матрицы А хранится в процессоре i . На первом шаге первая строка рассылается всем процессорам, после чего вычисления (*) могут выполняться параллельно процессорами Р 2 , … , P n . На втором шаге вторая строка приведенной матрицы рассылается из процессора P 2 процессорам P 3 , … , P n , и т. д. На k - ом шаге процессоры P 1 , … , P k -1 простаивают, P k выполняет нормировку и рассылку строки, P k +1 , … , P n — вычисляют по формуле (*). Рассылка строки выполняется за время ( n - k )( n - k ) T send , нормировка — за время ( n - k ) T mul . Параллельная модификация строк выполняется за время ( n - k ) T mul . Таким образом Т.е. при больших значениях n параллельная версия может работать медленнее последовательной версии. Это объясняется значительным объемом обмена данными между каждыми двумя шагами и уменьшением числа активных процессоров на 1 на каждом шаге. В ситуации, когда p n , проблема балансировки нагрузки в определенной степени смягчается. Будем считать, что n = mp , и строки 1, p +1, 2 p +1, … хранятся в процессоре 1, строки 2, p +2, 2 p +2, … — в процессоре 2, и т.д. Такой способ хранения называется циклической слоистой схемой. Для уменьшения простоя процессоров будем использовать опережающую рассылку: на k -м шаге ( k +1)-я строка рассылается сразу после того, как окончится ее модификация. Для обеспечения численной устойчивости алгоритма применим схему выбора главного элемента по строке. В этом случае получим При достаточно больших значениях n Альтернативой для хранения элементов по строкам является столбцовая циклическая схема хранения. В этом случае по процессорам распределяются не строки, а столбцы, что позволяет уменьшить объем передаваемых данных.

Главный элемент будем выбирать в активном столбце. На k -м шаге процессор, хранящий k - й столбец, выбирает главный элемент, вычисляет значения i > k ) и пересылает их на остальные процессоры. Далее каждый процессор вычисляет по формуле Передача данных происходит за время p ( n - k ) T send , вычисления — за время . Тогда При достаточно больших значениях n Зафиксируем p =4 и построим графики функции S p ( n ) при различных значениях отношения T send / T mul . Вычислим оптимальное число процессоров p ’ и построим график функции S p ( n ) при n =5000. Значение p ’ найдем из уравнения p ³ 1 Отсюда получаем Результаты эксперимента

Рассмотрим только решение нижнетреугольной системы (решение верхнетреугольной системы организуется аналогичным образом). В качестве алгоритма решения выберем алгоритм скалярных произведений. Его запись на псевдокоде выглядит следующим образом для i = 1 до n для j = 1 до i -1 b i = b i – L ij y j y i = b i / L ii Распараллеливание Пусть вычислительная система состоит из p процессоров, а размерность матрицы L равна n . Разобьем матрицу L на блоки в соответствии с рисунком Все блоки (кроме, может быть, двух последних) состоят из p блоков разбиваются вектора y и b . Далее, пусть блок S 1 расположен в памяти процессора P 1 , блоки L 2 и S 2 — в памяти процессора P 2 и т.д. Перед началом вычислений блок b k пересылается процессору P k , k = 1.. p . Рассмотрим первый шаг алгоритма решения системы. 1. Процессор P 1 решает систему S 1 y 1 = b 1 2. Найденный вектор y 1 пересылается процессорам P 2 , …, P p 3. Процессор P k ( k =2.. p ) пересчитывает вектор b k по формуле Остальные шаги алгоритма организуются аналогично.

Оценки производительности Будем считать, что матрица L уже распределена по процессорам. В противном случае выгоднее решать систему последовательным способом. На одном процессоре задача размерности n решается за время Оценим время выполнения k -ой итерации в системе из p процессоров Тогда Зафиксируем p =4 и построим графики функции S p ( n ) при различных значениях отношения T send / T mul . Вычислим оптимальное число процессоров p ’ и построим график функции S p ( n ) при n =1500. Значение p ’ найдем из уравнения p ³ 1 Отсюда получаем Результаты эксперимента К сожалению, для данной задачи не удалось получить приемлимых результатов. Дело в том, что на момент проведения экспериментов в сети было недостаточное количество компьютеров, поэтому параллельный алгоритм выполнялся даже дольше, чем последовательный. 2.3. QR -разложение Постановка задачи Рассмотрим ортогональное приведение матрицы A A = QR где Q —ортогональная, а R — верхнетреугольная матрицы. В качестве алгоритма решения задачи выберем преобразование Хаусхолдера . На последовательных компьютерах преобразование Хаусхолдера выполняется примерно в два раза медленне , чем LU -разложение.

Поэтому данный метод редко используется для решения СЛУ, однако он широко используется при вычислении собственных значений. Пусть A — матрица n n . Рассмотрим запись алгоритма на псевдокоде для k = 1 до n -1 , u T k = (0, …, 0, a kk – s k , a k +1, k , …, a nk ) a kk = s k для j = k +1 до n a j = g k u k a j , a j = a j - a j u k Распараллеливание Первый шаг алгоритма выглядит следующим образом 1. А всем процессорам 2. s 1 в процессоре 1 и переслать его всем процессорам. 3. А во всех процессорах в соответствии с алгоритмом Хаусхолдера . Остальные шаги организуются аналогичным образом.

Оценки производительности Пусть А — матрица размерности n n . И пусть система состоит из p процессоров, причем n = pm . Распределим матрицу А по процессорам в соответствии со столбцовой циклической схемой хранения.

Последовательный алгоритм выполняется за время Пересылка начальных данных выполняется за время n 2 T send , пересылка результата — за время n 2 T send . Поэтому Зафиксируем p =4 и построим графики функции S p ( n ) при различных значениях отношения T send / T mul . Построим график функции S p ( n ) при n =2000. Результаты эксперимента

Приведем две стандартные теоремы сходимости.

Теорема. Если матрица А имеет строгое диагональное преобладание или является неприводимой с диагональным преобладанием, то итерации метода Якоби сходятся при любом начальном приближении х 0 . Теорема. Если матрица A = D – B симметрична и положительно определена, то итерации метода Якоби сходятся при любом начальном приближении х 0 тогда и только тогда, когда матрица D + B положительно определена. Для проверки сходимости итерационных приближений будем использовать следующий критерий: где . Распараллеливание Пусть система состоит из р процессоров и матрицы разбиты на блоки следующим образом , Припишем блоки H 1 , d 1 , x k 1 процессору Р 1 , блоки H 2 , d 2 , x k 2 — процессору Р 2 и т.д. Тогда соответствующие компоненты вектора х можно вычислять параллельно по схеме: Шаг 1. Переслать на i - й процессор H i , d i , x 0 . Шаг 2. На к-й итерации i - й процессор выполняет следующие действия: · х i к+1 = H i x i k + d i · x i k+1 остальным процессорам · x k +1 · Шаг 3. Если хотябы один флаг не установлен, то перейти к следующей итерации (Шаг 2). Шаг 4. Передать результат.

Оценки производительности Для простоты будем считать, что n = pm . Оценим время выполнения одной итерации метода Якоби. Далее, пусть метод сошелся за k шагов. Тогда Из последней формулы несложно установить, что если число итераций , то каким бы большим небыло значение n , нам все равно не удастся достичь ускорения S p ³ 2. Зафиксируем p =4, k =500 и построим график функции S p ( n ) при различных значениях отношения T send / T mul . Кривая, которой соответствует значение T send / T mul =1000 никогда не поднимется выше 1, т.к. число итераций k слишком мало.