Прикладной специалист в процессе решения своей прикладной задачи, скорее всего, хотел бы знать как можно меньше об устройстве конкретной ЭВМ (но при этом полностью использовать ее возможности); не иметь проблем при переходе с одной машины на другую; пользоваться при программировании привычными средствами и не вникать в “программистсткие” тонкости, типа особенностей компиляторов, языков программирования, организации системных библиотек, вопросов эффективного их использования и т.п.

При разработке прикладных программ для параллельных ЭВМ такую комфортную модель обычно реализовать не удается. Основные причины этого - необходимость знать и использовать понятия из новой предметной области - области параллельных вычислений. Это и архитектура многопроцессорных систем (массивно-параллельные системы (MPP), симметричные мультипроцессорные системы (SMP), кластерные системы и т.п., используемые микропроцессоры, коммуникационная сеть, организация памяти - распределенная или общая); и модели параллельности (типы параллельных процессов и способы их задания, передача сообщений); и представление моделей параллельности в языковых средствах программирования (стандарты MPI или PVM для модели передачи сообщений или OpenMP для SMP-систем в модели параллельности с общей памятью).

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

Прежде чем начать обсуждение системы Норма и результатов ее применения, приведем классификацию типов сеток, используемых при решении задач математической физики (сайт Northeast Parallel Architectures Center, Syracuse University, www.npac.syr.edu/PROJECTS/PUB/nikos/ParGG.html):

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

В Норма-программе указываются только те правила (ограничения), которым должны удовлетворять значения переменных, при этом в записи отсутствуют понятие памяти и другие мотивы, связанные с архитектурой целевого компьютера, а также обычные элементы программ (управление, циклы, переприсваивание и т.п.). Язык Норма является языком с однократным присваиванием (single assignment language) - каждая переменная в программе может принимать значение только один раз, а индексы переменных имеют вид i ± D , D - целая константа. Эти ограничения (и некоторые другие) связаны с необходимостью решать в процессе трансляции задачу синтеза выходной программы, которая в более общих постановках может не иметь решения.

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

Язык Норма является специализированным языком - первоначально он предназначался для описания решения вычислительных задач на статических структурных сетках, однако область его применения оказалась значительно шире. В разделе 3 будет приведен пример применения языка Норма для решения задачи на структурной статической сетке, а в разделе 4 изложен опыт применения Нормы для решения задач на статических неструктурных сетках и статических вложенных сетках.

В настоящее время реализованы версии транслятора с языка Норма, которые позволяют получать выходные программы для распределенных параллельных вычислительных систем (модель памяти NORMA, параллельные вычисления с передачей сообщений), и параллельных вычислительных систем с общей памятью (модель памяти NUMA, параллельные вычисления над общей памятью) и последовательных ЭВМ. Выходные программы могут быть построены на языках Фортран-77, Фортран PVM, Фортран MPI, Фортран GNS [3], Фортран Convex [4].

Рассмотрим фрагмент программы решения системы нестационарных уравнений газовой динамики с одной пространственной переменной методом Годунова [5]. Расчетная схема для рассматриваемой задачи:

r _iⁿ = r _i^n-1 - (D t / h) [(RU) _i+1 ^n-1 - (RU) _i ^n-1],

r u _iⁿ = (r u) _i^n-1 - (D t / h) [(P+RU) _i+1 ^n-1 - (P+RU²) _i+1 ^n-1],

r e _iⁿ = (r e) _i^n-1 - (D t / h) [(RUE+PU) _i+1 ^n-1 - (RUE+PU) _i ^n-1],

p = r (g - 1) (e - u² / 2)

Здесь p - давление, u - скорость, r - плотность, e – полная энергия единицы массы газа. Из первого уравнения определяется r , из второго u, из третьего e, а давление досчитывается из уравнения состояния идеального газа, g - показатель адиабаты газа. Рассчитанный таким образом набор газодинамических параметров принимается за исходный для следующего шага по времени и повторяется цикл расчета для момента времени t _n+1 = t _n + D t.

Программа на Норме, соответствующая этой вычислительной схеме.

! Начало главной программы TEST203.

!Описание переменных: NTEND - число шагов по времени, DT - шаг по времени,

! HX - шаг по пространству, GAMMA - показатель адиабаты

! Раздел GOD21.

!Вспомогательные сетки и переменные

! Условие завершения итерации

Файл исходных данных norma.dat:

Чтобы использовать при расчетах уже имеющиеся данные или выводить значения величин во внешнюю среду, применяются описания соответственно входных (INPUT) или выходных (OUTPUT) величин. Такое описание означает, что подлежат вводу (выводу) значения всех величин, указанных в списке ввода (вывода). Порядок, в котором будут вводиться (выводиться) величины, при этом не задается - он определяется в процессе трансляции автоматически (что требует, конечно, специальной обработки транслятором входных файлов).

FOR <область> ASSUME <соотношение>.

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

При описании итерационного процесса вычислений (конструкция ITERATION … END ITERATION) в заголовке итерации указан индекс итерации NT и перечислены итерируемые величины, являющиеся результатом этого фрагмента вычислений. Начальные значения итерируемых переменных заданы в конструкции INITIAL … END INITIAL. Далее описан собственно итерационный процесс: сначала описано вычисление "распадных" величин PX, RX, UX, EX на сетках GRIDIN, GRIDL, GRIDR по правилу RIEMANN, а затем искомых величин R, U, E, P. Правило RIEMANN в данном фрагменте не конкретизировано: оно может быть определено в другой части Норма-программы, при помощи Фортран-программы (как это было сделано в рассматриваемом случае) или библиотечной программы. Условие завершения итерации описано в конструкции EXIT WHEN.

1. Применение системы Норма для решения задач на статических неструктурных и статических вложенных сетках на параллельных ЭВМ с раcпределенной памятью

Рассмотрим применение языка Норма для решения на распределенных системах конкретной задачи газовой динамики, для которой в [6] на основе аппроксимации кинетического уравнения переноса схемой с направленными разностями на ячейке с произвольным числом граней построена кинетически-согласованная разностная схема (КСРС) и приведена дискретизация Навье-Стоксовских потоков для расчета газодинамических течений на треугольных сетках. Полученные схемы применены для расчета внешних задач обтекания. Расчет итеративный, схема расчета явная: все исходные данные для расчета текущего слоя итерации берутся с предыдущего слоя. Таким образом, можно вычислять расчетные переменные в каждом узле сетки независимо от значений в других узлах, то есть можно говорить, что данная задача имеет уровень параллельности порядка n, где n - число узлов сетки.

При разработке Норма-программы для данной задачи возникают две основные проблемы: 1) как задать в языке Норма вычислительную формулу при зависимости одной величины от другой в случае использовании неструктурной сетки; 2) как распределить расчетные величины между процессорами таким образом, чтобы добиться сбалансированности вычислений и минимизировать время на обмены между процессорами. Рассмотрим эти проблемы более подробно.

Первая проблема определяется использованием нерегулярных сеток и неформально может быть сформулирована следующим образом. Пусть необходимо вычислить значение расчетной величины X в точке (i,j). При вычислениях на регулярных сетках аргументы для вычисления величин обычно находятся в точках с координатами вида: (i+k,j+l), где k,l -целочисленные константы. При вычислениях на нерегулярных сетках точки сетки нумеруются от 1 до N. Пусть точки, соседние к i-й точке сетки, имеют номера i1,i2,...,ik. Тогда общая структура сетки определяется матрицей neigh, состоящей из N строк, в k-й строке указываются номера точек, соседних к k-й точке. Например, если i-я строка содержит номера i1,i2,...,ik, и для вычисления значения величины X в точке i необходимо использовать значения X1 в соседних точках, то запись в программе выглядит неформально следующим образом (F - некоторая функциональная зависимость):

Запись таких выражений в явном виде в языке Норма не предусмотрена (отсутствует понятие массива и косвенной адресации).

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

r _i^k = r _i^k-1 - (dt / v _i ) (S _j - D _j + … ) (1),

Здесь W (i) - граница - определяет множество точек, соседних к точке i.

На языке Фортран такую запись обычно можно представить в следующем виде (vsn(i) - число соседей у точки с номером i):

В случае, когда такая программа должна исполняться на параллельной ЭВМ с общей памятью, достаточно сделать указание о том, что внешний цикл по i можно выполнять параллельно. Однако, необходимо учитывать, что переменная sum разделяется всеми витками цикла и поэтому необходимо указать, что она является локальной для каждого витка цикла. Например, в Фортране-Convex [4] это можно сделать следующим образом:

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

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

Один из способов преодоления этой трудности состоит в том, чтобы преобразовать исходную запись (1) к следующему виду:

r _i^k = r _i^k-1 - (dt / v _i ) (S1 _i,j - D1 _i,j + … ) (2),

В этом случае вместо значений S(j) и D(j) в каждой точке мы вводим вспомогательные массивы S1(i,j) и D1(i,j), которые строятся по следующему правилу: S1(i,j)=S(neigh(i,j)), то есть "размножаем" используемые значения.

Если выражение (1) мы не могли записать на Норме, то в последнем случае такая запись возможна и имеет вид:

FOR oi ASSUME ro =ro(it-1)-dt/v*SUM((oj) S1-D1),

где oi : (i=1..N), oj : (j=1..8)/j<vsn(i) - используемые области изменения индексов. При этом величины S1 и D1 описываются на области oij : (oi;oj):

VARIABLE S1,D1 DEFINED ON oij.

Здесь важно отметить следующее обстоятельство. Если при программировании записи (1) вычисленные значения величин S и D привязываются к конкретным точкам, то в случае (2) необходимо после завершения витка итерационного цикла предусмотреть распределение вычисленных значений по соответствующим местам в матрицах S1 и D1. Для этого и используются функции, выполнение которых обеспечивает обмен данными между процессорами. Здесь важно отметить также то, что на каждом шаге итерации обмен происходит между одними и теми же процессорами, так как сетка статическая - конфигурация сетки не изменяется на протяжении счета задачи. Конечно, нетрудно заметить, что при таком подходе существенно увеличивается размер используемой памяти. Однако, теперь для каждой точки известно, что аргументы, используемые для ее вычисления, находятся в том же процессоре, где находится и сама точка.

Вторая проблема связана с распределением расчетных точек по процессорам распределенной системы таким образом, чтобы минимизировать время на передачу сообщений между ними. Эта проблема в последнее время интенсивно изучается (например, [7-9]) и является самостоятельной задачей, не зависящей от средств и методов разработки программ. В общем случае время, затрачиваемое на обмены информацией между процесорами, определяется двумя характеристиками. Во-первых, это количество вызовов процедур обмена, то есть, если точки, помещенные в процессор, имеют соседство с точками, размещенными в других k процессорах, то тогда необходимо обращаться к процедуре обмена k раз. Второй характеристикой является объем передаваемой информации между процессорами.

Время, необходимое для выполнения операции отправки сообщения, в общем случае, имеет вид T=t*L+b,

где L - длина сообщения в байтах, t - чистое время передачи одного байта и b - накладные расходы на вызов процедуры передачи данных (latency time). Известно, что b>> t. Время T определяет накладные расходы, вызванные тем, что задача решается на распределенной вычислительной системе.

Минимизировать T можно как снижением числа обращений к процедурам обмена (уменьшение k), так и уменьшением общего количества передаваемой информации.

В системах [7-9] основное внимание уделяется сокращению общего объема передаваемой информации, причем, зачастую за счет увеличения числа обращений к процедурам обмена. К преимуществам таких систем следует отнести тот факт, что они работают очень быстро. Кроме того, уже появляются параллельные версии таких систем, которые позволяют использовать их для динамически изменяемых сеток.

В нашем случае мы использовали собственную программу, в которой основное внимание уделяется сокращению числа обращений к процедурам обмена. Например, для сетки, состоящей из 75900 точек удалось так распределить точки по процессорам, что каждый процессор, в процессе счета, обменивался данными только со своими двумя соседями. При этом общий объем передаваемой информации превосходил в два раза объем, полученный при разбиении с помощью системы METIS (версия 3.0). При использовании системы METIS число обращений к процедурам обмена варьировалось от 2 до7 для различных процессоров.

Ниже приведены времена счета на ЭВМ МВС-100, МВС-1000К, Parsytec-CC12 (12 процессоров), Parsytec-CC32 (24 процессора) для задачи моделирования дозвукового течения вязкого газа, рассматриваемой выше (зависимости от числа процессоров, ЭВМ и числа параллельных процессов).

Проводилось трехмерное численное моделирование развития крупномасштабной конвективной неустойчивости внутри протонейтронной звезды, возникшей в результате действия неравновесного процесса нейтронизации вещества. Система уравнений адиабатической гидродинамики с учетом гравитации решается с помощью явной консервативной TVD разностной схемы Годуновского типа со вторым порядком по пространству и времени. Для повышения точности результатов используется метод вложенных сеток. Расчетная область покрывается набором вложенных друг в друга сеток с одинаковым числом ячеек, но последовательно уменьшающимся в 2 раза абсолютным размером. Граница сетки n-уровня всегда совпадает с границей сетки n-1 уровня. Для выполнения условия устойчивости необходимо, чтобы двум шагам вычисления на сетке n-уровня соответствовал один шаг вычислений по времени на сетке n-1-го уровня. Пересчет всех величин в каждой ячейке на сетке n-1 уровня совершался простым усреднением 8 значений этих величин на сетке n - го уровня. Граничные значения, необходимые для вычисления величин на сетке n -го уровня получались в результате монотонной интерполяции по 27 значениям этих величин на сетке n-1 уровня. В вычислениях использовалось три уровня вложенных сеток (G0,G1,G2), каждая размерностью 56х56х56, с постоянным шагом по пространству ( схема сетки приведена ниже).

Сетка 3-го уровня с наименьшим шагом покрывает полностью начальное возмущение энтропии. Это позволяет на два порядка увеличить число расчетных ячеек в этой области, по сравнению со случаем одной сетки. Один шаг на сетке 0-го уровня соответствует 2 шагам на сетке 1-го уровня и четырем шагам на сетке 3-го уровня: всего получается 7 временных шагов.

На каждом временном шаге сначала определяются граничные значения для сетки Gi по значениям в соответствующих узлах сетки G(i-1). Далее вычисляются значения расчетных величин во внутренних узлах сетки Gi. В заключение полученные результаты усредняются и возвращаются в соответствующие узлы сетки G(i-1). Поскольку на сетке каждого уровня используется одна и та же гидродинамическая часть программы (в дальнейшем будем называть ее Hydra), то увеличение числа уровней приводит к большому росту общего времени расчета, что приводит к необходимости разработки параллельной программы.

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

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

С другой стороны, такое распределение ячеек сетки по процессорам приводит к трудностям, связанным с организацией вычислений в граничных ячейках и с возвратом усредненных значений на сетку более высокого уровня. Это связвано с тем, что в обоих случаях ячейки, для которых проводится расчет, и аргументы, используемые при этом, находятся, вобще говоря, в различных процессорах. Таким образом, для решения данной задачи надо запрограммировать вспомогательные алгоритмы, которые обеспечивали бы обмен необходимыми данными. Сетка i-го уровня накладывается на сетку уровня i-1, начиная с ячейки с номером (Mx,My,Mz), то есть в ячейке (Mx,My,Mz) сетки Gi-1 размещаются 8 ячеек сетки Gi. Значения Mx,My,Mz задаются в качестве параметров исходной задачи. Это позволяет до начала основного расчета определить, значения из каких ячеек необходимо передавать из одного процессора в другой, то есть для каждого процессора построить вспомогательные массивы, в которых будет храниться информация о том, в какой процессор (из какого процессора) и значения из каких ячеек надо передавать (принимать). В дальнейшем для обменов можно использовать уже готовую информацию из этих массивов, не проводя лишних вычислений.

Обратим внимание на вопрос, где надо проводить соответствующие вычисления. Рассмотрим решение этого вопроса на примере вычисления значений в граничных ячейках сетки G2 с использованием значений в ячейках сетки G1. Количество граничных ячеек приблизительно равно 12*N*N, число используемых аргументов - 4.25*N*N. Очевидно, что в этом случае вычисление значений в граничных ячейках надо проводить в процессорах, в которых эти ячейки находятся. Аналогично рассматривается вопрос о том, где вычислять усредненные значения. Очевидно, что количество усредненных значений в 8 раз меньше числа используемых аргументов. В этом случае усреднение проводится в процессорах, в которых находятся аргументы, а затем результаты рассылаются по необходимым процессорам.

1. И.Б.Задыхайло, К.Н.Ефимкин. Содержательные обозначения и языки нового поколения. Информационные технологии и вычислительные системы, N2, 1996, c. 46-58.
2. А.Н.Андрианов, А.Б.Бугеря, К.Н.Ефимкин, И.Б.Задыхайло. Норма. Описание языка. Рабочий стандарт. Препринт ИПМ им.М.В.Келдыша РАН, N120, 1995, 50 с.
3. В.А.Абрамова и др. Система программирования GNS. Препринт ИПМ им.М.В.Келдыша РАН, N71, 1997, 23с.
4. Convex Fortran User’s Guide. Order N DSW-038, Convex Press, Richardson, Texas, USA, October 1994, pp. 1-277.
5. А.Н.Андрианов, С.Б.Базаров, К.Н.Ефимкин. Решение задач газовой динамики методом Годунова на параллельных ЭВМ с помощью языка Норма. Препринт ИПМ им.М.В.Келдыша РАН, N 121, 1995, 22с.
6. И.В.Абалакин, А.В.Жохова, Б.Н.Четверушкин. Кинетически-согласованные разностные схемы на нерегулярных сетках. Математическое моделирование, т.9, N7, 1997, c.44-53.
7. L.Oliker. PLUM: Parallel Load Balancing for Adaptive Unstructered Meshes. Journal of Parallel and Distributed Computing 52, 150-177(1998).
8. G.Karypis, K.Schloegel, V.Kumar. ParMETIS - Parallel Graph Partitioning and Sparse Matrix Ordering Library. Version 2.0. 1998 ( http: //www.cs.umn.edu/ ~karypis).
9. C.Walshaw, M.Cross. Parallel Optimisation Algorithms for Multilevel Mesh Partitioning. Tech. Rep. 99/IM/44, Univ. Greenwich, London SE18 6PF, UK, February 1999.

старший научный сотрудник ИПМ им. М.В.Келдыша РАН, к.ф.-м.н.

e-mail: and@a5.kiam.ru

зав. отделом ИПМ им. М.В.Келдыша РАН, к.ф.-м.н.

e-mail: efi@a5.kiam.ru