95

Как написать программу на языке Норма

А.Н. Андрианов, К.Н. Ефимкин, С.Д. Устюгов

ИПМ им. М.В. Келдыша РАН

Содержание

• 1. Введение
• 2. Математическая постановка задачи об отражении наклонной ударной волны и метод ее решения
• 3. Основные конструкции языка Норма
• 3.1 Способ задания областей
• 3.2 Способ задания величин
• 3.3 Задание вычисления переменной
• 3.3.1 Оператор ASSUME
• 3.3.2 Оператор COMPUTE
• 3.3.3 Оператор COMPUTE в операторе ASSUME
• 3.3.4 Интерфейс с программами, написанными на Фортране
• 3.4 Итерационный цикл
• 4. Программа для распределенной системы: общие проблемы и их автоматическое решение в системе Норма
• 5. Разработка программы: от Фортрана к Норме
• 5.1 Общая схема вычислений исходной программы
• 5.2 Главная программа
• 5.3 Распределенная программа
• 5.3.1 Вычисление временного шага dt для итерации по времени
• 5.3.2 Вычисление потоков F и G
• 5.3.2.1 Вычисление величины eq
• 5.3.2.2 Вычисление величины qsr
• 5.3.2.3 Вычисление величин slp и csr
• 5.3.2.4 Вычисление величины ggf
• 5.3.2.5 Вычисление величины as, gmx, ed, f
• 5.3.3 Программа раздела mhdx для вычисления потока F
• 5.3.4 Завершение итерации по методу Рунге-Кутта
• 5.3.5 Вычисление величины qu на границах
• 6. Трансляция и запуск Норма-программы
• 7. Литература
• 8. Приложение 1. Текст последовательной Фортран-программы
• 9. Приложение 2. Текст Норма-программы
• 10. Приложение 3. Текст Фортран-программ, используемых в Норма-программе

1. Введение

Цель данного материала - достаточно подробно описать основные приемы разработки параллельной программы на языке Норма [1-3]. Ориентирован этот материал в первую очередь на прикладного специалиста, который сталкивается с необходимостью разработки параллельной программы для решения своей расчетной задачи.

Надо отметить, что в настоящее время имеется много различных методов и средств разработки параллельных программ (с ними можно познакомиться, например, в [4]). Эти способы требуют, в различной форме, достаточно сложной и трудоемкой работы программиста по выявлению параллелизма и представлению этого параллелизма в программе. Надежды на автоматическое распараллеливание последовательной программы пока не оправдываются: для параллельных вычислительных систем с общей памятью распараллеливающие компиляторы не обеспечивают эффективное распараллеливание, а для параллельных вычислительных систем с распределенной памятью проводятся лишь исследовательские работы в данном направлении.

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

Программирование на языке Норма также требует усилий - не должно создаваться иллюзии, что применение этого языка снимает все проблемы параллельного программирования, конечно, это не так. Самое очевидное соображение против использования этого языка - необходимость изучать этот новый язык. Однако усилия по написанию Норма-программы, по нашему мнению, оправданы - вся работа с программой на многопроцессорных вычислительных системах (изменение размеров сетки, изменение числа процессоров, изменение выходного языка) поддерживается системой, а проблемы по синтезу этой параллельной программы остаются только у Норма-компилятора, а не у программиста. Кроме того, имеется положительный опыт применения этого языка для решения практических задач математической физики (например, [5-9]).

В качестве исходной информации для разработки Норма-программы взята последовательная программа на Фортране 77 для численного решения задачи об отражении сильной ударной волны от стенки, предоставленная С.Д. Устюговым (ИПМ РАН). Следует отметить, что такой подход является не вполне правильным - разработку Норма-программы следует вести по расчетным формулам, а не по готовой последовательной программе, при этом система Норма позволяет получать по одной и той же Норма-программе как параллельную, так и последовательную реализацию. Перевод последовательной программы на язык Норма является трудоемким процессом -- требуется, по крайней мере, детальный анализ исходной последовательной программы, чтобы понять, а что же необходимо вычислить. При этом приходится "распутывать" информационные зависимости, часто не связанные с существом расчета, а возникшие как результат определенного стиля программирования и богатства возможностей Фортрана (экономия памяти, косвенная адресация, массовое использование COMMON-блоков и т.п.). Тем не менее, мы опишем процесс перехода от последовательной Фортран-программы к Норма-программе, поскольку эта ситуация является типичной - прикладной специалист чаще всего предпочитает именно такую форму представления необходимых расчетов.

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

2. Математическая постановка задачи об отражении наклонной ударной волны и метод ее решения

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

Начальные значения физических величин в расчётной области определялись с помощью соотношений Ренкина-Гюгонио для ударной волны с числом Маха при показателе адиабаты и имели вид

перед фронтом ударной волны и



за фронтом ударной волны.

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

где

, ,

и

- вектор консервативных переменных,

,

- вектора потоков,

- плотность,

- компоненты скорости,

- давление,

- полная энергия,

- удельная внутренняя энергия.

Эта система решалась с помощью явной консервативной конечноразностной TVD схемы Годуновского типа второго порядка по пространству и времени [10]

где - правый собственный вектор матрицы , - правый собственный вектор матрицы . Элементы вектора находятся по формуле

.

Значения есть собственные значения матрицы , соответствующие характеристичеким скоростям и определялись для некоторого среднего состояния по известным значениям и по методу Рое [11]. Функция есть энтропийная коррекция к и имеет вид [12]:

,

где - константа, равная в расчёте .

В вычислениях использовалась функция-лимитер

.

Численный поток определялся по аналогичным формулам.

Продвижение во времени от к осуществлялось методом Рунге-Кутта третьего порядка сохраняющего TVD свойство исходной конечноразностной схемы [13].

Шаг по времени определялся из условия Куранта

В вычислениях число Куранта бралось равным 0.5. Расчёты проводились на равномерной сетке размерностью .

3. Основные конструкции языка Норма

Специализированный декларативный язык нового поколения Норма [1-3] (разработка ИПМ им.М.В.Келдыша РАН) позволяет описывать решение широкого класса задач математической физики и при этом использовать привычные для прикладного специалиста понятия. Такое описание не ориентировано на конкретную архитектуру компьютера, поэтому оно предоставляет большие возможности для выявления естественного параллелизма и организации вычислений. Программа на языке Норма фактически является представлением математических расчетных формул, полученных прикладным специалистом.

Способ задания областей

Понятие области введено в языке Норма для представления понятия индексного пространства. Неформально, область - это множество значений, которые могут принимать (или принимают) индексы величин, используемых в расчетных формулах. Точнее, область - это совокупность целочисленных наборов {i₁,....,i_n }, n>0, i_j>0, j=1,...,n, каждый из которых задает координаты точки n-мерного индексного пространства. С каждым направлением (осью координат) n-мерного пространства задачи связывается уникальное имя - имя индекса (имя оси координат индексного пространства).

Следует отметить, что область определяет значения координат точек индексного пространства, а не значения расчетных величин в этих точках. Например, если требуется вычислить значения величины Y_i,j, i, j =1,...,n на некоторой сетке X_i,j, i,j =1,...n, заданной, например, формулой X_i,j =F(h,i,j), (F - заданная функция, h - заданный параметр), то следует:

описать область, состоящую из точек (i,j), i,j = 1,...n;

(2) описать на этой области величины X и Y;

(3) задать на этой области правило вычисления значений сетки X_i,j : X_i,j = F(h,i,j) и правило вычисления значений Y_i,j : Y_i,j = G(X_i,j) (функция G также некоторым способом задана).

В языке Норма область может иметь имя. Индексы областей специально не описываются - они вводятся при определении областей. Области используются при описании величин, определенных на области, при задании вычислений в операторах ASSUME, при описании входных или выходных величин, при задании областей фактических параметров в вызовах разделов или функций, в функциях редукции.

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

oi:(i=1..Nx). oj:(j=MyJ..Ny+1). om:(m=1..4).

Имя-одномерной-области может использоваться для ссылки на эту область. Имя-индекса есть индексная переменная, множество значений которой определяется диапазоном. Границами диапазона являются целые положительные константные выражения, построенные из целых констант, параметров области и арифметических операций. Требуется, чтобы в программе всем параметрам области были присвоены конкретные значения в описании параметров области. Для этого в языке предусмотрена конструкция DOMAIN PARAMETER, например, DOMAIN PARAMETER Nx=388, MyJ=4, Ny=122. Таким образом, области oi и oj представляют собой наборы точек 1..388 и 4..123 соответственно.

Многомерная область строится при помощи операции “;” произведения прямоугольных областей. Например, описания трехмерной области oijm и двумерной oijstep, заданных с помощью операции произведения одномерных областей, могут быть описаны следующим образом:

oijm:(oi;oj;om). oijstep:((i=3..Nx-2);(j=3..Ny-2)).

Способ задания величин

Скалярные величины (скаляры) и величины на области относятся к арифметическим величинам. Описание ставит в соответствие каждой арифметической величине уникальное в текущей программной единице имя величины, а также задает тип величины: REAL, INTEGER или DOUBLE (по умолчанию - тип REAL). Пример описания скаляров:

VARIABLE dt0,dk REAL.

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

Например, описание

VARIABLE qu DEFINED ON oijm

говорит о том, что определена величина qu на области oijm, то есть эта величина имеет индексы i,j,m и может принимать значения во всех точках области oijm.

Если индексное выражение у величины не содержит смещения, то оно может быть опущено; например, записи qu[i, j, m-1] и qu[m-1] эквивалентны.

Если значение некоторого индекса задается константой, необходимо явно указать, к какому направлению относится константа, например, обращение к величине qu_i,j,m с индексами i=5,j=3,m=2 выглядит следующим образом: qu[i=5,j=3,m=5].

Если необходимо значения одного индекса связать со значениями другого, необходимо сделать явное указание, например, диагональные элементы в плоскости m=2 можно определить следующим образом: qu[i,j=i,m=2].

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

Задание вычисления переменной

В Норме определены три вида операторов: скалярный оператор, оператор ASSUME и оператор COMPUTE. Операторы предназначаются для описания вычислительных действий, необходимых для решения задачи.

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

3.1.1 Оператор ASSUME

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

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

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

FOR ((i=3..Nx-2);(j=3..Ny-2);(m=1..4)) ASSUME qu=1.0.

во всех точках области, указанной после слова FOR, величине qu будет присвоено значение 1.0 (индексы у величины qu в записи соотношения опущены).

Если на одной и той же области надо вычислить несколько величин, можно задать несколько соотношений, например,

FOR onedomain ASSUME qu=1.0; e=p+v; u=w[i-1]/2.

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

og1:(i=3..Nx-2);(j=Ny-1..Ny);(m=1..4)).

og2:(i=3..Nx-2);(j=1..2);(m=1..4)).

Тогда оператор

FOR og1,og2 ASSUME qu=qu1.

описывает то же вычисление, что и два оператора

FOR og1 ASSUME qu=qu1.

FOR og2 ASSUME qu=qu1.

Язык Норма является языком с однократным присваиванием (single assignment language), величины могут принимать значения только один раз, переприсваивание значений невозможно по определению. Таким образом, некоторое внешнее сходство оператора ASSUME с оператором присваивания не должно вводить в заблуждение.

3.1.2 Оператор COMPUTE

Оператор COMPUTE вызова раздела в определенном смысле является аналогом понятия вызова подпрограммы в традиционных языках программирования и, по существу, является обобщением понятия соотношения, используемого в операторе ASSUME, так как дает возможность получать несколько величин-результатов.

При задании оператора COMPUTE надо задать значения величин, которые являются исходными данными для вычисления (фактические in-параметры) и имена величин, которые являются результатами вычисления (фактические out-параметры). Например, оператор

COMPUTE step(gam1,gamma,c0,dx,qu ON oijm RESULT umin ON oijstep).

является запросом на вычисление с именем step (подробности которого описаны в разделе с именем step), с in-параметрами gam1,gamma,c0,dx,qu (последний параметр - все значения величины qu на области oijm), по которым будут вычислены значения величины-результата umin (величина umin будет вычислена во всех точках области oijstep).

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

Фактическими out-параметрами могут быть скаляры, величины с индексами (быть может, заданными правилом умолчания индексов), величины на областях.

Вычисление раздела возможно, если все фактические in-параметры уже вычислены (приняли некоторые значения), побочный эффект в языке Норма невозможен.

3.1.3 Оператор COMPUTE в операторе ASSUME

Оператор COMPUTE может использоваться для задания сложного соотношения в операторе ASSUME, например

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,qu ON oijm/(i=i,j=j) RESULT umin).

Эта запись означает, что в каждой точке области oijstep надо провести вычисление step1, (подробности которого описаны в разделе с именем step1), с in-параметрами gam1,gamma,c0,dx,qu (последний параметр - значения величины qu на динамической области), по которым будет вычислено значение единственного результата - величины umin.

Фактический in-параметр qu ON oijm/(i=i,j=j) задает значения величины qu на динамической области, то есть области, которая меняется при различных значениях индексов из заголовка oijstep оператора ASSUME. Неформально, для каждого фиксированного значения индексов (i,j) области oijstep, например, i=i0, j=j0, при вызове раздела step1 в качестве фактического параметра передается рабочая величина R(1:4), такая, что R(1)=qu(i0,j0,1); R(2)=qu(i0,j0,2); R(3)=qu(i0,j0,3); R(4)=qu(i0,j0,4).

В целом, осуществление вычисления step1 в каждой точке области oijstep позволяет вычислить значения величины-результата umin во всех точках области oijstep, то есть вычисление step из предыдущего пункта и вычисление step1 в операторе ASSUME дадут одинаковый результат - величина umin будет вычислена во всех точках области oijstep. Разделы step и step1 при этом, конечно, будут отличаться.

Таким образом, в п.3.3.2 и п.3.3.3 приведены два способа описания вычислений с помощью оператора вызова раздела.

Отметим, что если вызов раздела находится в теле оператора ASSUME, то параметрами-результатами не могут быть скаляры и величины на статических областях (которые не меняются при различных значениях индексов из заголовка оператора ASSUME) - это заведомо приводит к переприсваиванию, что недопустимо в языке Норма.

3.1.4 Интерфейс с программами, написанными на Фортране

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

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

Итерационный цикл

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

ITERATION <список итерируемых величин> ON <индекс итерации>.

BOUNDARY

<операторы присваивания граничных значений итерируемым величинам>

END BOUNDARY

INITIAL <индекс итерации>=0:

<операторы присваивания начальных значений итерируемым величинам>

END INITIAL

<операторы>

EXIT WHEN (<условие выхода из итерации>).

END ITERATION <индекс итерации>.

Граничные значения итерируемой величины задаются при помощи обычных операторов языка Норма внутри блока BOUNDARY...END BOUNDARY. Эти значения считаются неизменными на протяжении всей итерации и определены на каждом шаге итерации.

Начальные значения для итерируемых величин задаются операторами языка Норма при помощи блока INITIAL...END INITIAL.

Индекс итерации может использоваться в списке индексов, задаваемом для итерируемой переменной, то есть фактически можно считать, что итерируемая величина в пределах итерации имеет дополнительный итерационный индекс по фиктивному направлению. Правило умолчания индекса справедливо и для итерационного индекса. Например, запись qu[itime-1], где qu - итерируемая величина, itime - индекс итерации, задает значения величины qu с предыдущего шага итерации. Обычно в пределах итерации сохраняются значения итерируемых переменных для двух шагов итерационного цикла: текущего и предыдущего.

4. Программа для распределенной системы: общие проблемы и их автоматическое решение в системе Норма

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

При построении фрагмента параллельной распределенной программы, реализующей параллельное вычисление на области O некоторой переменной X, определенной на области D, необходимо определить:

1) номера процессоров, которые должны выполнять данный фрагмент программы,

2) границы изменения параметров циклов, обеспечивающих обход области O,

3) описание массива X,

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

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

Проиллюстрируем на примере, в чем заключается решение задач 1-4). Пусть надо вычислить следующие соотношения:

1. U(i,j)=i+j,для i=1..N; j=1..N.

2. W(i,j)=U(i,j)-1, для i=1..N; j=1..N.

3. V(i,j)= U(i-1,j)+1, для i=25..N-25; j=1..N.

Пусть N=100 и распараллеливание следует проводить на линейку из 10 процессоров по направлению i. На языке Норма эти вычисления описываются следующим образом:

! Описание областей Oi, Oi1, Oj

Oi:(i=1..N). Oi1:(i=25..N-25). Oj:(j=1..N).

! Область Oij есть декартово произведение (i=1..N)X(j=1..N)

Oij:(Oi;Oj).

! Область Oij1 есть декартово произведение (i=25..N-25)X(j=1..N)

Oij1:(Oi1;Oj).

! Определение параметра областей

DOMAIN PARAMETERS N=100.

! Описание переменных на областях

VARIABLE U,V,W DEFINED ON Oij.

! Число процессоров по направлению i равно 10.

DISTRIBUTION INDEX i=1..10, j=1.

! Расчетные формулы (индексы без смещений можно опускать)

FOR Oij ASSUME U=i+j; W=U-1.

FOR Oij1 ASSUME V=U[i-1]+1.

Фрагмент выходной программы, который будет построен компилятором с языка Норма при трансляции на язык Фортран-77 с использованием библиотеки MPI, представлен ниже.

C Вычисление параметров циклов в зависимости от номера процессора

IN01 = 1

IK01 = 10

IN02 = 1

IK02 = 10

IF(IPRWWW.EQ.3) THEN IN02=5

IF(IPRWWW.EQ.8) THEN IK02=5

DO j=1,100

DO i=IN01,IK01

C Выражение i+j в процессоре с номером IPR

U(i,j)=((IPRWWW-1)*10+i)+j

ENDOO

ENDDO

C Требуются вычисленные значения из процессоров с номерами 3..7

IF(IPRWWW.LT.3.OR.IPRWWW.GT.7.OR.JPRWWW.NE.1) GOTO 4

DO j=1,100

RM1(j)=U(10,j)

ENDDO

CALL MPI_SEND(RM1,100,MPI_REAL, ITASKJ(JPRWWW+(IPRWWW+1-1)),

>2,MPI_COMM_WORLD,IER1)

4 CONTINUE

DO j=1,100

DO i=IN01,IK01

W(i,j)=U(i,j)-1

ENDOO

ENDDO

C Используются вычисленные значения в процессорах с номерами 4..8

IF(IPRWWW.LT.4.OR.IPRWWW.GT.8.OR.JPRWWW.NE.1) GOTO 7

CALL MPI_RECV(RM2,100,MPI_REAL,ITASKJ(JPRWWW+(IPRWWW-2)),

>2,MPI_COMM_WORLD,STTS1,IER1)

DO j=1,100

U(0,j)=RM2(j)

ENDDO

7CONTINUE

IF(IPRWWW.LT.3.OR.IPRWWW.GT.8) GOTO 8

DO j=1,100

DO i=IN02,IK02

V(i,j)= U(i-1,j)+1

ENDOO

ENDDO

CONTINUE

Видно, что эта программа менее проста и понятна, чем программа на Норме, однако это скорее всего не должно беспокоить пользователя - данная программа построена автоматически компилятором. Поясним, используя текст этой программы, как решены задачи 1-4). Идентификатор IPRWWW обозначает номер процессора из линейки процессоров и принимает значения от 1 до 10, идентификатор JPRWWW=1 (представление линейки процессоров).

Определение номеров процессоров. Пусть некоторое вычисление необходимо выполнить при значениях индекса i, принимающего значения в диапазоне от до . В языке Норма значения и известны на этапе трансляции, поэтому компилятор по значению () и числу точек, распределенных на процессор, определяет номер процессора () в котором находится точка (). Таким образом, диапазон процессоров, участвующих в этом вычислении, задается сегментом [, ]. Процессоры, номера которых находятся за пределами этого сегмента, не участвуют в реализации данного вычисления. В нашем примере операторы, реализующие соотношения 1-2, выполняются во всех 10 процессорах. Соотношение 3 вычисляется при значениях индекса i в диапазоне от 25 до 75. Точка i=25 принадлежит процессору с номером 3, а точка i=75 принадлежит процессору с номером 8. Поэтому оператор, реализующий соотношение 3, выполняется только в процессорах, номера которых лежат в диапазоне [3,8] - это условие реализуется последним условным оператором фрагмента программы.

Определение границ изменения параметров циклов. Аналогично тому, как определяется диапазон процессоров, реализующих заданное соотношение (вычисление), определяются начальные и конечные значения параметра цикла i, при которых надо выполнить оператор, реализующий соотношение. В нашем случае, операторы 1 и 2 выполняются (по индексу i) от IN01 до IK01, а оператор 3 - от IN02 до IK02. При этом во всех 10 процессорах, значения IN01=1, IK01 = 10. Более сложный вид имеют значения IN02 и IK02. В процессоре с номером 3 значение IN02 = 5, во всех остальных IN02 = 1. В процессоре с номером 8 значение IK02 = 5, а в остальных IK02 = 10.

Описания массивов. Построение описаний для величин, определенных на областях, проводится достаточно простым образом. В последовательной программе переменным U,V,W соответствуют описания

DIMENSION U(100,100), V(100,100), W(100,100)

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

DISTRIBUTION INDEX i=1..10, j=1

то есть "разрезаны" на 10 процессоров по частям, состоящим из 10100 элементов. Кроме этого, поскольку в формуле для вычисления V требуются значения U[i-1], необходимо выделить память ("теневую" грань) для приема значений U[i-1] из соседнего левого процессора. Это делается за счет расширения влево диапазона при описании массива U. Таким образом, описания для распределенной программы переменных U,V,W выглядят так:

DIMENSION U(0:10,100), V(10,100), W(10,100)

Генерация операторов обмена. Интерфейс передачи сообщений MPI предоставляет разнообразные средства реализации взаимодействия в модели передачи сообщений. В рассматриваемой программе взаимодействие между процессорами реализовано при помощи подпрограмм MPI_SEND (отправить сообщение) и MPI_RECV (принять сообщение), которые имеют следующие параметры:

CALL MPI_SEND(<адрес посылки>,<длина посылки>,<тип посылаемых элементов>,

<имя процесса-получателя>,<тэг>,<коммуникатор >,<код возврата>)

CALL MPI_RECV (<адрес приема>,<длина приема>, <тип принимаемых элементов>,

<имя процесса-отправителя>,<тэг>, <коммуникатор>,<статус>,<код возврата>.

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

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

ITASKJ(n,m) - матрица имен параллельных процессов, участвующих в вычислении (ее размеры n и m фактически совпадают с параметрами распределения, заданными в конструкции DISTRIBUTION INDEX, поэтому в данном примере описание матрицы имен имеет вид INTEGER ITASKJ(10,1);

IPRWWW, JPRWWW - значения индексов матрицы ITASKJ, определяющие имя каждого из процессов (другими словами, каждый процесс "знает", что его имя - ITASKJ(IPRWWW,JPRWWW)).

Таким образом, вызов

CALL MPI_SEND (RM1,100,MPI_REAL, ITASKJ(JPRWWW+IPRWWW),

2,MPI_COMM_WORLD,IER1)

приводит к пересылке 100 элементов массива RM1 типа REAL соседу справа. Посылку осуществляет каждый из процессов с номерами JPRWWW=1, IPRWWW=3..7.

Вызов

CALL MPI_RECV (RM2,100,MPI_REAL,ITASKJ(JPRWWW+IPRWWW-2),

2,MPI_COMM_WORLD,STTS1,IER1)

приводит к приему 100 элементов типа REAL в массив RM2 от соседа слева. Прием осуществляет каждый из процессов с номерами JPRWWW=1, IPRWWW=4..8. Тэг со значением 2 позволяет “связать” прием и отправку сообщения (это важно, когда один и тот же процессор посылает одному и тому же процессору несколько сообщений - разные тэги позволяют однозначно идентифицировать эти сообщения).

Определение параметров операторов MPI_SEND и MPI_RECV проводится автоматически на этапе компиляции. Место в выходной программе, где должны стоять эти операторы, определяется также компилятором из следующих соображений. Оператор MPI_SEND ставится непосредственно за оператором, в котором вычисляется требуемое для передачи значение. Оператор MPI_RECV ставится непосредственно перед оператором, в котором используется принимаемое (требуемое) значение. Это позволяет надеяться на перекрытие времени счета и обмена данными. В нашем примере, сразу после вычисления значений переменной U стоит оператор записи значений величины U в рабочий массив RM1 и оператор отправки этого массива MPI_SEND. При этом, как отмечено выше, отправляют значения не все процессоры. Соответствующий оператор приема MPI_RECV расположен перед оператором, где используется значение U[i-1,j]. Соответствие обеспечивается наличием уникального тэга 2, причем генерацию тэгов также поддерживает компилятор. Отметим, что компилятор с языка Норма поддерживает несколько методов расстановки операторов обмена сообщениями.

В заключение несколько слов о схеме распараллеливания вычислений при трансляции Норма-программы.

Пусть распределенная вычислительная система представлена как матрица процессорных элементов. Как и при задании матриц можно считать, что каждый элемент такой вычислительной системы характеризуется своими координатами: номером строки (i) и номером столбца (j).

В программе на языке Норма программист сам определяет, точки какого индексного направления будут распределяться по процессорам системы, и задает число процессоров, на которых будет решаться исходная задача. Направление "разрезания" и число процессоров задается, например, описанием следующего вида: DISTRIBUTION INDEX i=1..6;j=1..5. Такая запись означает, что все области, в описании которых используются направления i или j, будут распределяться на 6 процессоров по направлению i и на 5 процессоров по направлению j.

Для эффективного распараллеливания при задании индексных направлений i, j в конструкции DISTRIBUTION INDEX обычно используются следующие два критерия:

1) максимальное число точек по индексному направлению (что приводит к распределению между процессорами наиболее "тяжелой работы"),

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

Разрезание исходных данных проводится с учетом выбранных направлений. Обозначим верхнюю границу по направлению i для величины через . Пусть . Используемая процедура разрезания основывается на размещении постоянного числа компонентов всех массивов в каждом процессоре (по фиксированному направлению). Пусть по направлению i используется m процессоров, а по направлению j - n процессоров. Тогда число компонентов массива по направлению i (соответственно j) задается значением

().

Пусть координата процессора по направлению i есть iprwww, по направлению j суть jprwww. Тогда в процессоре с номером (iprwww, jprwww) вычисляются значения по направлению i от (iprwww-1)ipoints+1 до iprwwwipoints, а по направлению j от (jprwww-1)jpoints+1 до jprwwwjpoints. В нашем примере процессоры образуют матрицу размерностью 65. Исходная матрица значений имеет размерность . Таким образом, каждый процессор проводит расчет величины qu размерностью

.

5. Разработка программы: от Фортрана к Норме

5.1 Общая схема вычислений исходной программы

Текст исходной программы на Фортране приведен в Приложении 1. Вычисления в этой программе организованы по следующей схеме:

Присваивание начальных значений расчетным переменным

Итерационный цикл по времени

Вычисление временного шага для итерации по времени

Итерационный цикл Рунге-Кутта (3 шага)

Вычисление потока F

Вычисление потока G

Вычисление величины qu

Конец итерации Рунге-Кутта

Вычисление значений величины qu на границах

Конец итерации по времени

Запись вычисленных значений величины qu в файл

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

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

Главная программа

В исходной программе значения переменной qu многократно перевычисляются, так как Фортран допускает переприсваивание значений переменных. В программе на Норме в таких случаях будут вводиться дополнительные “копии” величин для исключения переприсваивания (в выходной программе компилятор оптимизирует использование требуемой из-за введения лишних переменных лишней памяти). Для определения начальных значений величины qu на первом этапе вводится величина qu0. Главная программа на Норме выглядит следующим образом:

MAIN PART GDST.

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

oijm:((i=1..Nx);(j=1..Ny);(m=1..4)). ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU DEFINED ON oijm.

VARIABLE x1 DEFINED ON ok1. VARIABLE x2 DEFINED ON ok2.

VARIABLE gamma,gam1,c0,dt0,an,fi,pi,dx REAL.

pi=3.1415926.

gamma=1.4.

gam1=gamma-1.0.

c0=0.5.

an=1.0/6.0.

fi=pi/3.0.

COMPUTE init(pi,gam1,fi,an

RESULT dx,x1 ON ok1, x2 ON ok2, QU0 ON oijm,dt0).

COMPUTE GD(pi,gamma,gam1,c0,fi,an,dx,dt0,x1 ON ok1,x2 ON ok2,QU0 ON oijm

RESULT QU ON oijm).

COMPUTE writer(x1 ON ok1, x2 ON ok2, QU ON oijm).

END PART.

В этой программе, кроме присваивания значений скалярным переменным, заданы вызовы трех вычислений - init, GD, и writer. Вычисление init является подпрограммой, написанной на Фортране и используемой в рамках интерфейса Нормы и Фортрана. Вычисление GD является разделом языка Норма и будет описано на Норме. Вычисление writer также является подпрограммой, написанной на Фортране. Рассмотрим эти три части программы.

Вызов раздела

COMPUTE init(pi,gam1,fi,an

RESULT dx,x1 ON ok1, x2 ON ok2, QU0 ON oijm,dt0).

с входными параметрами pi, gam1, fi приводит к вычислению значений величины qu0 во всех точках области oijm (i=1..Nx; j=1..Ny ;m=1..4), значений векторов x1 и x2; значения шага по пространству dx и начального значения dt0. Так как расчет начальных данных будет выполняться в главной программе и будет выполняться на одном процессоре, то можно просто использовать соответствующий код последовательной Фортран-программы (см. Приложение 1), оформив его в виде подпрограммы init, который практически совпадает с начальным фрагментом текста исходной программы SHATL на Фортране.

SUBROUTINE init(pi,gam1,fi,an,dx,x1,x2,QU,dt0)

INCLUDE 'TASK.PAR'

real QU(Nx,Ny,4)

real x1(Nx),x2(Ny)

db=1.0

dx=db/(Ny-1)

x1(1)=-2.*dx

DO I=2,Nx

x1(I)=x1(I-1)+dx

ENDDO

x2(1)=-2.*dx

DO I=2,Ny

x2(I)=x2(I-1)+dx

ENDDO

DO J=1,Ny

DO I=1,Nx

IF(x1(i).LT.AN+x2(J)/TAN(fi)) THEN

QU(I,J,1)=8.

QU(I,J,2)=8.25*SIN(PI/4.)*QU(I,J,1)

QU(I,J,3)=-8.25*COS(PI/4.)*QU(I,J,1)

QU(I,J,4)=116.5/gam1+(QU(I,J,2)**2+QU(I,J,3)**2)/2./QU(I,J,1)

ELSE

QU(I,J,1)=1.4

QU(I,J,2)=0.

QU(I,J,3)=0.

QU(I,J,4)=1./gam1

ENDIF

ENDDO

ENDDO

dt0=0.0

RETURN

END

В данной подпрограмме, а также в приводимых ниже подпрограммах, используется ссылка на файл 'TASK.PAR'. Этот файл содержит следующие описания:

PARAMETER(NX=388,NY=122,ipoints=65,jpoints=25)

Cipoints=[Nx/iproc], jpoints=[NY/jproc]

PARAMETER(iproc=6,jproc=5)

INTEGER ITASKJ(iproc*jproc)

COMMON /NORMAMESSPASS/ITASKJ,IPRWWW,JPRWWW,IVVVVJ,IWWWWJ

В описаниях задаются значения основных параметров задачи: Nx, Ny - размерности сеток по направлению i и по направлению j; iproc (jproc) - число процессоров по направлению i(j); ipoints(jpoints) - число точек сетки в одном процессоре по направлению i(j). Последняя строка содержит описание стандартного (для Норма-системы) COMMON-блока. В частности, iprwww указывает номер процессора по направлению i, а jprwww указывает номер процессора по направлению j, а массив ITASKJ задает все процессы, участвующие в параллельной работе. Неформально, этот файл обеспечивает интерфейс программ на Фортране, построенных компилятором с языка Норма, и программ на Фортране, написанных вручную.

Оператор

COMPUTE GD(pi,gamma,gam1,c0,fi,an,dx,dt0,x1 ON ok1,x2 ON ok2,QU0 ON oijm

RESULT QU ON oijm).

задает вызов раздела GD, который будет определен как распределенный (предназначенный для параллельного выполнения) и в котором будут описаны основные вычисления задачи. Исходными параметрами для раздела GD являются значения величины qu0, вектора x1 и x2 и ряд констант. Результатом вычисления раздела являются значения расчетной величины qu во всех точках области определения.

Оператор

COMPUTE writer(x1 ON ok1, x2 ON ok2, QU ON oijm).

задает вызов раздела writer для записи вычисленных значений величины qu в файл. Раздел writer, как и раздел init, практически совпадает с заключительным фрагментом текста исходной программы SHATL на Фортране и оформлен в виде подпрограммы.

SUBROUTINE writer(x1,x2,QU)

INCLUDE 'TASK.PAR'

real QU(Nx,Ny,4)

real x1(Nx)

real x2(Ny)

open(40,file='REZ')

WRITE(40,*) 'VARIABLES="X","Y","D","U","V"'

WRITE(40,*) 'ZONE T="XX",I=118,J=384,F=POINT'

DO I=3,NX-2

DO J=3,NY-2

TT=QU(I,J,1)

UU=QU(I,J,2)/TT

VV=QU(I,J,3)/TT

WRITE(40,*) X1(I),X2(J),TT,UU,VV

ENDDO

ENDDO

close(40)

RETURN

END

Распределенная программа

Теперь перейдем к определению распределенного раздела GD. Величина qu вычисляется в результате итерационного процесса (цикла) по времени itime, причем на каждом шаге итерации вычисление этой величины разбивается на 2 этапа. На первом этапе вычисление проводится с помощью 3-х шагов итерации по методу Рунге-Кутта (по переменной irk). На втором этапе перевычисляются значения на нижней и верхней границах исходной области.

Величина qu остается неизменной в областях (i=1..2);(j=1..Ny);(m=1..4) и (i=Nx-1..Nx);(j=1..Ny);(m=1..4) - эти области в программе на языке Норма обозначены соответственно og1 и og2. Поэтому значения величины qu в этих областях необходимо определять в конструкции граничных значений (BOUNDARY...ENDBOUNDARY) итерации по itime. В областях og3: ((i=3..Nx-3);(j=1..2);(m=1..4)) и og4: ((i=3..Nx-3);(j=Ny-1..Ny);(m=1..4)) величина qu перевычисляется на каждом шаге итерации по времени itime. При этом на каждой области вычисления проводятся по различным формулам и реализуются в разделах downbound и upbound.

	Ny					og4
	Ny-1
		og1			oin			og2
J
	2					og3
	1
		1	2						Nx-1	Nx

Рис.1. Схема расчетных областей задачи

В области oin: ((i=3..Nx-3);(j=3..Ny-2);(m=1..4)) вычисление величины qu проводиться с помощью 3-х шагов итерации методом Рунге-Кутта. В остальных точках области при проведении итерации методом Рунге-Кутта значения величины qu остаются неизменными, поэтому в программе на Норме эти значения также задаются в конструкции BOUNDARY...ENDBOUNDARY итерации по irk. Схема расчетных областей с точки зрения организации итерационных вычислений по методу Рунге-Кутта приведены на рис.2.

	Ny					og4
	Ny-1	Граница для итерации Рунге-Кутта
				Область итеративных вычислений Рунге-Кутта
		og1				oin				og2
J
	2					og3
	1
		1	2						Nx-1	Nx

Рис.2. Схема расчетных областей для итерации Рунге-Кутта

Программа распределенного раздела GD приведена ниже. Строки, начинающиеся с символа "!", в программе на языке Норма являются комментарием.

PART GD.

! входные параметры для распределенной программы

pi,gamma,gam1,c0,fi,an,dx,dt0,x1,x2,QU0

! результаты вычисления раздела

RESULT QU

BEGIN

DOMAIN PARAMETERS Nx=388,Ny=122.

DISTRIBUTION INDEX i=1..6,j=1..5.

oi:(i =1..Nx). oi3nm2:(i=3..Nx-2). oi2nm2:(i=2..Nx-2).

oj:(j =1..Ny). oj3nm2:(j=3..Ny-2). oj2nm2: (j =2..Ny-2).

m4:(m =1..4). m5: (m =1..5).

! i=1,Nx j=1,Ny,m=1,4

oijm:(oi;oj;m4).

!i=3,Nx-2 j=3,Ny-2,m=1,4 i=3,Nx-2 j=3,Ny-2

oin:(oi3nm2;oj3nm2;m4). oijstep:(oi3nm2;oj3nm2).

!i=3,Nx-2 j=1,Ny,m=1,4

oin1:(oi3nm2;oj;m4).

!i=1,2 j=1,Ny,m=1,4 i=Nx-1,Nx j=1,Ny; m=1,4

og1:((i=1..2);oj;m4). og2:((i=Nx-1..Nx);oj;m4).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og3:(oi3nm2;(j=1..2);m4). og4:(oi3nm2;(j=Ny-1..Ny);m4).

og33:(oi3nm2;(j=1..2)). og44:(oi3nm2;(j=Ny-1..Ny)).

!i=3,Nx-2 j=1,2;m=1,4 i=3,Nx-2, j=Ny-1,Ny; m=1,4

og5:(oi3nm2;(j=1..2);m4). og6:(oi3nm2;(j=Ny-1..Ny);m4).

!i=1,Nx j=1,Ny i=1,Nx j=1,Ny m=1,5

oprx:(oi;oj). oprw:(oprx;m5).

!i=2,Nx-2 j=3,Ny-2 m=1,4 i=3,Nx-2 j=2,Ny-2 m=1,4

oif:(oi2nm2;oj3nm2;m4). oig:(oi3nm2;oj2nm2;m4).

ok1:(k=1..Nx). ok2:(k=1..Ny).

VARIABLE QU0,QU,QU2 DEFINED ON oijm.

VARIABLE umin DEFINED ON oijstep.

VARIABLE x1 DEFINED ON ok1.

VARIABLE x2 DEFINED ON ok2.

VARIABLE f DEFINED ON oif.

VARIABLE g DEFINED ON oig.

VARIABLE prx DEFINED ON oprx.

VARIABLE prw DEFINED ON oprw.

VARIABLE pi,gam1,gamma,c0,dt,dt0,dt0it,dk,dx,dtr,an,fi REAL.

VARIABLE A1,A2,A9,B9,C9,D9 REAL.

A9=0.25. B9=0.75. C9= 0.33333333. D9=0.66666666. dk=0.15.

ITERATION QU,dt0it ON itime.

BOUNDARY

!i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

FOR og1;og2 ASSUME QU=QU0.

END BOUNDARY

INITIAL itime=0:

dt0it=dt0.

! i=3..Nx-2; j=1..Ny;m=1..4

FOR oin1 ASSUME QU=QU0.

END INITIAL

!i=3..Nx-2; j=3..Ny-2

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,QU[itime-1] ON oijm/(i=i,j=j) RESULT umin).

dt =MIN((oijstep)umin).

dt0it=dt0it[itime-1]+dt.

OUTPUT dt(FILE='tst').

dtr=dt/dx.

ITERATION QU2 ON irk.

BOUNDARY

! i=1..2; j=1..Ny;m=1..4 & i=Nx-1..Nx; j=1..Ny;m=1..4

! i=3..Nx-2; j=1..2;m=1..4 & i=3..Nx-2; j=Ny-1..Ny;m=1..4

FOR og1;og2;og3;og4 ASSUME QU2=QU[itime-1].

END BOUNDARY

INITIAL irk=0:

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME QU2=QU[itime-1].

END INITIAL

! i=1..Nx;j=1..Ny

FOR oprx ASSUME prx=gam1*(QU2[irk-1,m=4]

-(QU2[irk-1,m=2]**2+QU2[irk-1,m=3]**2)/2.0/QU2[irk-1,m=1]).

! i=1,Nx;j=1,Ny; m=1,3

FOR oprw/m=1..3 ASSUME prw=QU2[irk-1].

! i=1,Nx;j=1,Ny; m=4

FOR oprw/m=4 ASSUME prw=QU2[irk-1]+prx.

! i=1,Nx;j=1,Ny; m=5

FOR oprw/m=5 ASSUME prw=gamma*prx/QU2[irk-1,m=1].

! i=2..Nx-2;j=3..Ny-2;m=1..4

COMPUTE mhdx (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT f ON oif).

! i=3..Nx-2;j=2..Ny-2;m=1..4

COMPUTE mhdy (gam1,prw ON oprw, prx ON oprx, QU2[irk-1] ON oijm

RESULT g ON oig).

COMPUTE coeff (irk, A9,B9,C9,D9 RESULT A1,A2).

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

QU2=A1*QU[itime-1]+A2*(QU2[irk-1]-dtr*(f-f[i-1]+g-g[j-1])).

EXIT WHEN (irk=3).

END ITERATION irk.

! i=3..Nx-2; j=3..Ny-2;m=1..4

FOR oin ASSUME

COMPUTE absqu (QU2 ON oin/(i=i,j=j),i,j,an,x1 ON ok1

RESULT QU ON oin/(i=i,j=j)).

! i=3..Nx-2; j=1..2;m=1..4

FOR og33 ASSUME

COMPUTE downbound (QU ON og3/(i=i,j=3),i,x1 ON ok1,an,pi,gam1

RESULT QU ON og3/(i=i,j=j)).

! i=3..Nx-2; j= Ny-1,Ny;m=1..4

FOR og44 ASSUME

COMPUTE upbound (i,j,x1 ON ok1, x2 ON ok2,dt0it,fi,an,pi,gam1

RESULT QU ON og4/(i=i,j=j)).

EXIT WHEN (dt0it[itime]> dk).

END ITERATION itime.

END PART.

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

Величины prx и prw в исходной Фортан-программе определялись по одним и тем же формулам, как при вычислении потока f, так и при вычислении потока g. Разница состояла только в том, что вычисление проводилось при различных границах операторов цикла. В первом случае диапазон имел вид i=1,...,Nx; j=3,...,Ny-2 а во втором случае - i=3,...,Nx-2; j=1,Ny. В нашей программе мы вычисляем их один раз на области i=1,Nx;j=1,Ny.

Рассмотрим некоторые части этой программы более подробно. В начале программы приведены описания необходимых областей (эти описания определяются областями значений, которые принимают индексы в исходной программе на Фортране) и описания расчетных величин. Далее описываются необходимые вычисления.

5.1.1 Вычисление временного шага dt для итерации по времени

Значение временного шага dt зависит от значений величины qu в точках с координатами по оси i от 3 до Nx-2; а по оси j - от 3 до Ny-2.

Вычисление шага разбивается на две фазы. Cначала во всех точках рассматриваемой области вычисляются значения величины umin. Далее осуществляется поиск минимального значения среди элементов umin. Первая часть вычислений оформлена в виде подпрограммы step1. Для поиска минимального значения используется стандартная функция языка Норма (MIN). Все проблемы, связанные с вычислением этой функции на распределенной системе берет на себя транслятор с языка Норма.

В нашем случае вычисление величины umin выполняется в подпрограмме, в которой вычисляется один элемент величины umin. Вызов подпрограммы производится в операторе ASSUME на области oijstep:

FOR oijstep ASSUME

COMPUTE step1(gam1,gamma,c0,dx,QU[itime-1] ON oijm/(i=i,j=j) (*)

RESULT umin).

Выражение qu[itime-1] ON oijm/(i=i,j=j) означает, что в качестве параметра при вызове раздела step1, при фиксированных значениях параметров цикла - i=i0,j=j0, передается рабочий вектор R, состоящий из 4-х компонентов. При этом, значения соответствующих компонентов вектора R вычисляются по следующему правилу: R(1)=qu(i0,j0,1); R(2)=qu(i0,j0,2); R(3)=qu(i0,j0,3); R(4)=qu(i0,j0,4).

Возможен и другой подход к организации вычислений значений величины umin. При этом подходе в подпрограмме вычисляются все значения величины umin во всех точках области oijstep. Такое вычисление задается оператором:

COMPUTE step2(gam1,gamma,c0,dx,QU ON oijm RESULT umin ON oijstep). (**)

Ниже в таблице представлены описания двух подпрограмм, реализующих рассматриваемое вычисление величины umin.

SUBROUTINE step1(gam1,gamma,c0,dx,QU,umin) real QU(4) DTM=1.E+10 PRR=GAM1*(QU(4)-(QU(2)**2+ QU(3)**2)/2./QU(1)) CL=SQRT(GAMMA*PRR/QU(1)) UCSX=QU(2)/QU(1) UCSY=QU(3)/QU(1) DTMX=1./(ABS(UCSX)+CL) DTMY=1./(ABS(UCSY)+CL) umin=C0*DX*AMIN1(DTM,DTMX,DTMY) RETURN END	SUBROUTINE step2(gam1,gamma,c0,dx,QU,umin) INCLUDE 'TASK.PAR' real QU(ipoints,jpoints,4) real umin(ipoints,jpoints) in=1 ik=ipoints jn=1 jk=jpoints if(iprwww.eq.1)in=3 if(iprwww.eq.iproc)ik=NX-(iproc-1)*ipoints if(jprwww.eq.1)jn=3 if(jprwww.eq.jproc)jk=NY-(jproc-1)*jpoints DTM=1.E+10 do I=in,ik do J=jn,jk PRR=GAM1*(QU(I,J,4)-(QU(I,J,2)**2+ QU(I,J,3)**2)/2./QU(I,J,1)) CL=SQRT(GAMMA*PRR/QU(I,J,1)) UCSX=QU(I,J,2)/QU(I,J,1) UCSY=QU(I,J,3)/QU(I,J,1) DTMX=1./(ABS(UCSX)+CL) DTMY=1./(ABS(UCSY)+CL) umin(i,j)=C0*DX*AMIN1(DTM,DTMX,DTMY) enddo enddo RETURN END

Разница между подпрограммами step1 и step2 определяется следующими факторами.

Подпрограмма step1 вычисляет одно скалярное значение umin, не зависит от значений переменных по индексам распределения i и j, всю организацию распределенных вычислений берет на себя компилятор при трансляции оператора ASSUME (*).

Подпрограмма step2 имеет более сложный вид, так как она вычисляет значения umin на области, зависящей от индексов распределения i и j, зависит от значений переменной qu с индексами распределения i и j, а организация распределенных вычислений, сделанная компилятором при трансляции оператора (**) и отраженная в файле 'TASK.PAR', должна быть учтена при написании тела подпрограммы. Действительно, эта подпрограмма выполняется на каждом процессоре и поэтому границы используемых массивов по индексам распределения i и j зависят от размерности используемой сетки и числа процессоров, используемых для решения задачи (эти параметры содержатся в файле 'TASK.PAR'). Кроме того, существенным является вопрос определения границ для циклов, в которых вычисляются значения величины umin. В исходной программе диапазоны изменения параметров циклов имеет вид: i=3..Nx-2, j=3..Ny-2. При распределении на процессор по ipointsjpoints точек диапазоны изменения параметров циклов имеют вид i=1..ipoints, j= 1..jpoints для всех процессоров, кроме 4-х угловых (с координатами (1,1),(iproc,1),(1,jproc),(iproc,jproc) в матрице процессоров). Первые операторы подпрограммы step2 как раз и вычисляют эти значения - здесь фактически вручную программируются действия, которые для подпрограммы step1 будет сделана компилятором.