Файл: 12. Параллельные методы решения дифференциальных уравнений в частных производных.pdf

12. Параллельные методы решения дифференциальных
уравнений в частных производных
12. Параллельные методы решения дифференциальных уравнений в частных производных........................................................................................................................................... 1 12.1. Последовательные методы решения задачи Дирихле ................................................. 2 12.2. Организация параллельных вычислений для систем с общей памятью..................... 3 12.2.1. Использование OpenMP для организации параллелизма...................................... 3 12.2.2. Проблема синхронизации параллельных вычислений........................................... 4 12.2.3. Возможность неоднозначности вычислений в параллельных программах .......... 7 12.2.4. Проблема взаимоблокировки .................................................................................... 8 12.2.5. Исключение неоднозначности вычислений ............................................................. 9 12.2.6. Волновые схемы параллельных вычислений........................................................ 11 12.2.7. Балансировка вычислительной нагрузки процессоров......................................... 16 12.3. Организация параллельных вычислений для систем с распределенной памятью.. 17 12.3.1. Разделение данных.................................................................................................. 17 12.3.2. Обмен информацией между процессорами........................................................... 18 12.3.3. Коллективные операции обмена информацией .................................................... 20 12.3.4. Организация волны вычислений............................................................................. 21 12.3.5. Блочная схема разделения данных........................................................................ 22 12.3.6. Оценка трудоемкости операций передачи данных ............................................... 24 12.4. Краткий обзор раздела ................................................................................................... 25 12.5. Обзор литературы ........................................................................................................... 26 12.6. Контрольные вопросы..................................................................................................... 26 12.7. Задачи и упражнения ...................................................................................................... 27
Дифференциальные уравнения в частных производных представляют собой широко применяемый математический аппарат при разработке моделей в самых разных областях науки и техники. К сожалению, явное решение этих уравнений в аналитическом виде оказывается возможным только в частных простых случаях, и, как результат, возможность анализа математических моделей, построенных на основе дифференциальных уравнений, обеспечивается при помощи приближенных численных методов решения. Объем выполняемых при этом вычислений обычно является значительным и использование высокопроизводительных вычислительных систем является традиционным для данной области вычислительной математики. Проблематика численного решения дифференциальных уравнений в частных производных является областью интенсивных исследований (см., например, Березин и Жидков
(1966), Тихонов и Самарский (1977), Fox, et al. (1988)).
Рассмотрим в качестве учебного примера проблему численного решения задачи Дирихле для
уравнения Пуассона, определяемую как задачу нахождения функции
)
,
( y
x
u
u
=
, удовлетворяющей в области определения уравнению:
D
⎪
⎩
⎪
⎨
⎧
∈
=
∈
=
+
,
)
,
(
),
,
(
)
,
(
,
)
,
(
),
,
(
0 2
2 2
2
D
y
x
y
x
g
y
x
u
D
y
x
y
x
f
y
u
x
u
δ
δ
δ
δ
и принимающей значени
)
, y на гран е
0
я
x
g
иц
(
D об ас
D
( f и л ти
g являются функциями, задаваемыми при постановке задачи). Подобная модель может быть использована для описания установившегося течения жидкости, стационарных тепловых полей, процессов теплопередачи с внутренними источниками тепла и деформации упругих пластин и др. Данный пример часто используется в качестве учебно-практической задачи при изложении возможных способов организации эффективных параллельных вычислений (см. Немнюгин и Стесик (2002), Корнеев (1999),
Pfister (1998)).
Для простоты изложения материала в качестве области задания функции далее будет использоваться единичный квадрат:
D
)
,
( y
x
u
}
1
,
0
:
)
,
{(
≤
≤
∈
=
y
x
D
y
x
D

12.1. Последовательные методы решения задачи Дирихле
Одним из наиболее распространенных подходов численного решения дифференциальных уравнений является
метод конечных разностей (метод сеток) (см., например, Березин и Жидков (1966),
Тихонов и Самарский (1977), Pfister (1995)). Следуя этому подходу, область решения представляется в виде дискретного (как правило, равномерного) набора (
сетки) точек (узлов). Так, например, прямоугольная сетка в области может быть задана в виде (рис. 12.1):
D
D
⎩
⎨
⎧
+
=
+
≤
≤
=
=
=
),
1
/(
1
,
1
,
0
,
,
:
)
,
{(
N
h
N
j
i
jh
y
ih
x
y
x
D
i
i
j
i
h
где величина задает количество узлов по каждой из координат области
N
D
Обозначим оцениваемую при подобном дискретном представлении аппроксимацию функции в точках через
. Тогда, используя пятиточечный шаблон (см. рис. 12.1) для вычисления значений производных, мы можем представить уравнение Пуассона в конечно-разностной
форме:
)
,
( y
x
u
)
,
(
j
i
y
x
ij
u
ij
ij
j
i
j
i
j
i
j
i
f
h
u
u
u
u
u
=
−
+
+
+
+
−
+
−
2 1
,
1
,
,
1
,
1 4
Данное уравнение может быть разрешено относительно
:
ij
u
)
(
25 0
2 1
,
1
,
,
1
,
1
ij
j
i
j
i
j
i
j
i
ij
f
h
u
u
u
u
u
−
+
+
+
=
+
−
+
−
Разностное уравнение, записанное в подобной форме, позволяет определять значение по известным значениям функции в соседних узлах используемого шаблона. Данный результат служит основой для построения различных итерационных схем решения задачи Дирихле, в которых в начале вычислений формируется некоторое приближение для значений
, а затем эти значения последовательно уточняются в соответствии с приведенным соотношением. Так, например, метод
Гаусса-Зейделя для проведения итераций уточнения использует правило:
ij
u
)
,
( y
x
u
ij
u
)
(
25 0
2 1
1 1
,
1
,
,
1
,
1
ij
k
k
k
k
k
ij
f
h
u
u
u
u
u
j
i
j
i
j
i
j
i
−
+
+
+
=
−
−
+
−
+
−
, по которому очередное k-ое приближение значения вычисляется по последнему k-ому приближению значений и и предпоследнему (k-1)-ому приближению значений и
. Выполнение итераций обычно продолжается до тех пор, пока получаемые в результате итераций изменения значений не станут меньше некоторой заданной величины (требуемой точности вычислений). Сходимость описанной процедуры (получение решения с любой желаемой точностью) является предметом всестороннего математического анализа (см., например, Березин и Жидков (1966), Тихонов и Самарский
(1977), Pfister (1995)), здесь же отметим, что последовательность решений, получаемых методом сеток, равномерно сходится к решению задачи Дирихле, а погрешность решения имеет порядок
ij
u
j
i
u
,
1
−
1
,
−
j
i
u
j
i
u
,
1
+
1
,
+
j
i
u
ij
u
2
h
(i,j)
(i
+
1,j)
(i-1,j)
(i,j
+
1)
(i,j-1)
{
{
{
{
{
{
{
{
{
{
• • • • • • •
{
{
• • • • • • •
{
{
• • • • • • •
{
{
• • • • • • •
{
{
• • • • • • •
{
{
• • • • • • •
{
{
• • • • • • •
{
{
{
{
{
{
{
{
{
{
Рис. 12.1.
Прямоугольная сетка в области D (темные точки представляют внутренние узлы сетки, нумерация узлов в строках слева направо, а в столбцах - сверху вниз)
Рассмотренный алгоритм
(метод
Гаусса-Зейделя) на псевдокоде, приближенном к алгоритмическому языку С++, может быть представлен в виде:
// Алгоритм 12.1
// Метод Гаусса-Зейделя do { dmax = 0; // максимальное изменение значений u
2

3
for ( i=1; i }
} while ( dmax > eps );
Алгоритм 12.1.
Последовательный алгоритм Гаусса-Зейделя
(напомним, что значения при индексах
ij
u
1
,
0
,
+
= N
j
i
являются граничными, задаются при постановке задачи и не изменяются в ходе вычислений).
Рис. 12.2.
Вид функции в примере для задачи Дирихле
)
,
( y
x
u
Для примера на рис. 12.2 приведен вид функции
, полученной для задачи Дирихле при следующих граничных условиях:
)
,
( y
x
u
⎪
⎪
⎪
⎩
⎪⎪
⎪
⎨
⎧
=
+
−
=
+
=
−
=
−
∈
=
,
1
,
200 100
,
1
,
200 100
-
,
0
,
200 100
,
0
,
200 100
,
)
,
(
,
0
)
,
(
x
y
y
x
x
y
y
x
D
y
x
y
x
f
Общее количество итераций метода Гаусса-Зейделя составило 210 при точности решения
1 0
=
eps
и
(в качестве начального приближения величин использовались значения, сгенерированные датчиком случайных чисел из диапазона [-100,100]).
100
=
N
ij
u
12.2. Организация параллельных вычислений для систем с общей
памятью
Как следует из приведенного описания, сеточные методы характеризуются значительной вычислительной трудоемкостью:
2 1
kmN
T
=
, где есть количество узлов по каждой из координат области
,
- число операций, выполняемых методом для одного узла сетки, - количество итераций метода до выполнения условия остановки.
N
D
m
k
12.2.1. Использование OpenMP для организации параллелизма
Рассмотрим возможные способы организации параллельных вычислений для сеточных методов на многопроцессорных вычислительных системах с общей памятью. При изложении материала будем

предполагать, что имеющиеся в составе системы процессоры обладают равной производительностью, являются равноправными при доступе к общей памяти и время доступа к памяти является одинаковым
(при одновременном доступе нескольких процессоров к одному и тому же элементу памяти очередность и синхронизация доступа обеспечивается на аппаратном уровне). Многопроцессорные системы подобного типа обычно именуются симметричными мультипроцессорами (
symmetric multiprocessors,
SMP).
Обычный подход при организации вычислений для подобных систем – создание новых параллельных методов на основе обычных последовательных программ, в которых или автоматически компилятором, или непосредственно программистом выделяются участки не зависимых друг от друга вычислений. Возможности автоматического анализа программ для порождения параллельных вычислений достаточно ограничены, и второй подход является преобладающим. При этом для разработки параллельных программ могут применяться как новые алгоритмические языки, ориентированные на параллельное программирование, так и уже имеющиеся языки программирования, расширенные некоторым набором операторов для параллельных вычислений.
Оба указанных подхода приводят к необходимости значительной переработки существующего программного обеспечения, и это в значительной степени затрудняет широкое распространение параллельных вычислений. Как результат, в последнее время активно развивается еще один подход к разработке параллельных программ, когда указания программиста по организации параллельных вычислений добавляются в программу при помощи тех или иных внеязыковых средств языка программирования – например, в виде директив или комментариев, которые обрабатываются специальным препроцессором до начала компиляции программы. При этом исходный операторный текст программы остается неизменным, и по нему в случае отсутствия препроцессора компилятор может построить исходный последовательный программный код. Препроцессор же, будучи примененным, заменяет директивы параллелизма на некоторый дополнительный программный код (как правило, в виде обращений к процедурам какой-либо параллельной библиотеки).
Рассмотренный выше подход является основой
технологии OpenMP (см., например, Chandra, et al.
(2000)), наиболее широко применяемой в настоящее время для организации параллельных вычислений на многопроцессорных системах с общей памятью. В рамках данной технологии директивы параллелизма используются для выделения в программе
параллельных областей (parallel regions), в которых последовательный исполняемый код может быть разделен на несколько раздельных командных
потоков (threads). Далее эти потоки могут исполняться на разных процессорах вычислительной системы.
В результате такого подхода программа представляется в виде набора последовательных
(
однопотоковых) и параллельных (многопотоковых) участков программного кода (см. рис. 12.3).
Подобный принцип организации параллелизма получил наименование "
вилочного" (fork-join) или
пульсирующего параллелизма. Более полная информация по технологии OpenMP может быть получена в в разделе 5, в дополнительной литературе (см., например, Chandra, et al. (2000)) или в информационных ресурсах сети Интернет; в данном разделе возможности OpenMP будут излагаться в объеме, необходимом для демонстрации возможных способов разработки параллельных программ для рассматриваемого учебного примера решения задачи Дирихле.
12.2.2. Проблема синхронизации параллельных вычислений
Первый вариант параллельного алгоритма для метода сеток может быть получен, если разрешить произвольный порядок пересчета значений
. Программа для данного способа вычислений может быть представлена в следующем виде:
ij
u
//Алгоритм 12.2
// Параллельный метод Гаусса-Зейделя (первый вариант)
omp_lock_t dmax_lock;
omp_init_lock (dmax_lock);
do { dmax = 0; // максимальное изменение значений u
#pragma omp parallel for shared(u,n,dmax) \
private(i,temp,d)
for ( i=1; i#pragma omp parallel for shared(u,n,dmax) \
private(j,temp,d)
for ( j=1; j4

5
d = fabs(temp-u[i][j])
omp_set_lock(dmax_lock);
if ( dmax < d ) dmax = d;
omp_unset_lock(dmax_lock);
} // конец вложенной параллельной области
} // конец внешней параллельной области
} while ( dmax > eps );
Алгоритм 12.2.
Первый вариант параллельного алгоритма Гаусса-Зейделя
Следует отметить, что программа получена из исходного последовательного кода путем добавления директив и операторов обращения к функциям библиотеки OpenMP (эти дополнительные строки в программе выделены темным шрифтом, обратная наклонная черта в конце директив означает продолжение директив на следующих строках программы).
Рис. 12.3.
Параллельные области, создаваемые директивами OpenMP
Как следует из текста программы, параллельные области в данном примере задаются директивой
parallel for
, являются вложенными и включают в свой состав операторы цикла for. Компилятор, поддерживающий технологию OpenMP, разделяет выполнение итераций цикла между несколькими потоками программы, количество которых обычно совпадает с числом процессоров в вычислительной системе. Параметры директивы shared и private определяют доступность данных в потоках программы – переменные, описанные как shared, являются общими для потоков, для переменных с описанием private создаются отдельные копии для каждого потока, которые могут использоваться в потоках независимо друг от друга.
Наличие общих данных обеспечивает возможность взаимодействия потоков. В этом плане разделяемые переменные могут рассматриваться как
общие ресурсы потоков и, как результат, их использование должно выполняться с соблюдением
правил взаимоисключения (изменение каким-либо потоком значений общих переменных должно приводить к блокировке доступа к модифицируемым данным для всех остальных потоков). В данном примере таким разделяемым ресурсом является величина

900 358 20,61 358 295,03 0,07 358 11,03 1,87 1000 351 25,59 351 366,16 0,07 351 13,69 1,87 2000 367 106,75 367 1585,84 0,07 367 56,63 1,89 3000 370 243,00 370 3598,53 0,07 370 128,66 1,89
(k – количество итераций, T – время в сек., S – ускорение)
Оценим полученный результат. Разработанный параллельный алгоритм является корректным, т.е. обеспечивающим решение поставленной задачи. Использованный при разработке подход обеспечивает достижение практически максимально возможного параллелизма – для выполнения программы может быть задействовано вплоть до процессоров. Тем не менее, результат не может быть признан удовлетворительным – программа будет работать медленно и ускорение вычислений от использования нескольких процессоров окажется не столь существенным (результаты экспериментов говорят скорее всего о замедлении вычислений). Основная причина такого положения дел – чрезмерно высокая
синхронизация параллельных участков программы. В нашем примере каждый параллельный поток после усреднения значений должен проверить (и возможно, изменить) значение величины dmax.
Разрешение на использование переменной может получить только один поток – все остальные
2
N
ij
u
Вр ем я вы по лн ения
Параллельное выполнение
Ожидание
Последовательное выполнение
Процессоры (потоки)
1 2 3 4 5 6 7 8
Рис. 12.4.
Пример возможной схемы выполнения параллельных потоков при наличии синхронизации (взаимоисключения) потоки должны быть блокированы. После освобождения общей переменной управление может получить следующий поток и т.д. В результате необходимости синхронизации доступа многопотоковая параллельная программа превращается фактически в последовательно выполняемый код, причем менее эффективный, чем исходный последовательный вариант, т.к. организация синхронизации приводит к дополнительным вычислительным затратам – см. рис. 12.4. Следует обратить внимание, что, несмотря на идеальное распределение вычислительной нагрузки между процессорами, для приведенного на рис. 12.4 соотношения параллельных и последовательных вычислений, в каждый текущий момент времени (после момента первой синхронизации) только не более двух процессоров одновременно выполняют действия, связанные с решением задачи. Подобный эффект вырождения параллелизма из-за интенсивной синхронизации параллельных участков программы обычно именуются
сериализацией (serialization).
Как показывают выполненные рассуждения, путь для достижения эффективности параллельных вычислений лежит в уменьшении необходимых моментов синхронизации параллельных участков программы. Так, в нашем примере мы можем ограничиться распараллеливанием только одного внешнего цикла for. Кроме того, для снижения количества возможных блокировок применим для оценки максимальной погрешности многоуровневую схему расчета: пусть параллельно выполняемый поток первоначально формирует локальную оценку погрешности dm только для своих обрабатываемых данных (одной или нескольких строк сетки), затем при завершении вычислений поток сравнивает свою оценку dm с общей оценкой погрешности dmax.
Новый вариант программы решения задачи Дирихле имеет вид:
//Алгоритм 12.3
6