НЫСШЕЕ
ОБРАЗОВАНИЕ
А, А, Малявко
ПАРАЛЛЕЛЬНОЕ ПРОГРАММИРОВАНИЕ НА ОСНОВЕ ТЕХНОЛОГИЙ OPENMP, CUDA, OPENCL, MPI
3-е издание
УМО во
РЕКОМЕНДУЕТ
Еюрэйт
HiAUTiriLiTID
А. А. Малявко
ПАРАЛЛЕЛЬНОЕ ПРОГРАММИРОВАНИЕ НА ОСНОВЕ ТЕХНОЛОГИЙ OPENMP,CUDA,OPENCL, MPI
УЧЕБНОЕ ПОСОБИЕ ДЛЯ ВУЗОВ 3-е издание, исправленное и дополненное
Рекомендовано Учебно-методическим отделом высшего образования в качестве учебного пособия для студентов высших учебных заведений, обучающихся по инженерно-техническим направлениям
Книга доступна на образовательной платформе «Юрайт» urait.ru, а также в мобильном приложении «Юрайт.Библиотека»
Москва • Юрайт • 2022
https//и га it. ги
УДК 004.2(075.8)
ББК32.973я73
М21
Автор:
Малявко Александр Антонович — доцент, кандидат технических наук, доцент кафедры вычислительной техники факультета автоматики и вычислительной техники Новосибирского государственного технического университета.
Рецензенты:
Гунько А. В. — доцент, кандидат технических наук, доцент кафедры автоматики факультета автоматики и вычислительной техники Новосибирского государственного технического университета;
Ильиных С. П. — доцент, кандидат технических наук, доцент кафедры вычислительной техники факультета автоматики и вычислительной техники Новосибирского государственного технического университета.
Малявко, А. А.
М21 Параллельное программирование на основе технологий орешпр, cuda, opencl, mpi: учебное пособие для вузов / А. А. Малявко. — 3-е изд., испр. и доп. — Москва: Издательство Юрайт, 2022. — 135 с. — (Высшее образование). — Текст: непосредственный.
ISBN 978-5-534-14116-0
Курс содержит введение в технологии параллельного программирования для вычислительных систем различных классов. В него включены краткие сведения об архитектурах параллельных систем, рассмотрены понятия и виды параллелизма, описаны технологии параллельного программирования ОрепМР для систем с общей памятью, CUDA иОрепСБ для графических процесссоров и гетерогенных компьютеров, а также интерфейс передачи сообщений MPI для систем с распределенной памятью.
Соответствует актуальным требованиям федерального государственного образовательного стандарта высшего образования.
Для студентов высших учебных заведений, обучающихся по инженерно-техническим направлениям. Может использоваться студентами, магистрантами и аспирантами других специальностей при изучении родственных дисциплин, а также преподавателями смежных дисциплин
УДК 004.2(075.8)
ББК32.973я73 Все права защищены. Никакая часть данной книги не может быть воспроизведена в какой бы то ни было форме без письменного разрешения владельцев авторских прав.
ISBN 978-5-534-14116-0	© Малявко А. А., 2017
© Малявко А. А., 2021, с изменениями
© ООО «Издательство Юрайт», 2022
httpsI/и га it. ги
Оглавление
Предисловие...........................................5
Тема 1. Архитектура параллельных вычислительных систем................................................7
Тема 2. Введение в параллельное программирование......12
2.1.	Распараллеливание циклов.....................14
2.2.	Распараллеливание гнезд циклов...............19
2.3.	Преобразования гнезд циклов..................21
2.4.	Проблемы разработки параллельных программ....24
Тема 3. ОрепМР: параллельные вычисления на системах с общей памятью..........................29
3.1.	Модель параллельной программы ОрепМР.........29
3.2.	Директивы управления регионами и потоками....33
3.3.	Директивы распределения вычислений...........35
3.4.	Директивы синхронизации......................41
3.5.	Переменные окружения.........................43
3.6.	Библиотека функций ОрепМР....................44
Контрольные вопросы и задания......................46
Тема 4. CUDA: параллельные вычисления на графическом процессоре............................48
4.1.	Архитектура CUDA.............................48
4.2.	Аппаратные компоненты CUDA...................49
4.3.	Программные компоненты CUDA..................54
4.4.	CUDA-расширение языка C/C++..................56
4.5.	Технология разработки CUDA-программ..........60
4.6.	Библиотека CUDA-runtime, основные функции....66
Контрольные вопросы и задания.....................68
Тема 5. OpenCL: параллельные вычисления на гетерогенном компьютере...........................70
5.1.	Модель платформы.............................70
5.2.	Модель памяти OpenCL.........................71
3
5.3.	Модель исполнения OpenCL-программ..............73
5.4.	Модель структуры OpenCL-программ...............75
5.5.	Пример OpenCL-программы........................77
Контрольные вопросы и задания.......................83
Тема 6. MPI: параллельные вычисления на системах с распределенной памятью...............................85
6.1.	Основные понятия интерфейса MPI................86
6.2.	Функции инициализации и определения окружения..89
6.3.	Обмен сообщениями типа «точка-точка»...........90
6.4.	Коллективные операции взаимодействия процессов.96
6.5.	Производные типы и упаковка/распаковка данных..103
6.6.	Управление группами ветвей и коммуникаторами..108
6.7.	Управление виртуальными топологиями...........113
6.8.	Удаленный доступ к памяти (односторонние взаимодействия)....................................117
6.9.	Использование библиотеки МРЕ для анализа процессов взаимодействия ветвей программы....................125
Контрольные вопросы и задания......................131
Рекомендуемая литература..............................133
Предисловие
В курсе содержится краткое введение в теорию и практику параллельного программирования для параллельных систем различных классов (с общей памятью, с распределенной памятью, графические процессоры) и основные сведения о наиболее популярных технологиях поддержки параллельных вычислений.
Первая тема включает в себя краткие сведения об архитектурах многопроцессорных систем, необходимые для понимания существа параллельных вычислений.
Во второй теме рассматриваются: понятие параллелизма, основные виды параллелизма, описывается общая техника распараллеливания программ, а также специфические проблемы, сопровождающие разработку и отладку параллельных программ.
Третья тема посвящена описанию технологии параллельного программирования ОрепМР для систем с общей памятью.
В четвертой теме рассматривается архитектура CUDA как совокупность программно-аппаратных средств, обеспечивающих неграфические параллельные вычисления на графических процессорах производства компании NVIDIA.
Пятая тема посвящена изучению технологии OpenCL параллельного программирования для гетерогенных, т. е. содержащих разнотипные компоненты, вычислительных систем.
Шестая тема содержит подробное описание интерфейса передачи сообщений MPI, являющегося в настоящее время общепризнанным стандартом для программирования параллельных вычислительных систем с распределенной памятью.
Материал тем 3—6, содержащих описание конкретных технологий параллельного программирования, сопровождается контрольными вопросами.
В результате изучения дисциплины студент должен:
знать
•	компоненты программно-технических архитектур параллельных вычислительных систем;
5
•	виды параллелизма, уровни распараллеливания;
•	модель параллельной программы в технологии ОрепМР;
•	структуру памяти графического процессора и временные характеристики ее составных частей;
•	модели платформы, памяти, исполнения и структуры параллельной программы гетерогенного компьютера;
•	модель параллельной программы для вычислительной системы с распределенной памятью;
уметь
•	выявлять информационные зависимости между итерациями циклических участков программы;
•	применять директивы ОрепМР для создания параллельных регионов и распределения вычислительной работы между потоками;
•	применять рекомендуемую технологией CUDA стратегию параллельного программирования графического процессора;
•	применять рекомендуемую технологией OpenCL стратегию параллельного программирования гетерогенного компьтера;
•	использовать основные группы функций интерфейса MPI для инициации библиотеки, получения сведений о ее окружении;
владеть
•	способами преобразования циклов для ликвидации информационных зависимостей между итерациями;
•	способами борьбы с гонками данных, блокировками и тупиковыми ситуациями;
•	навыками разработки, компиляции и отладки параллельных программ с использованием технологии ОрепМР;
•	навыками разработки, компиляции и отладки параллельных программ с использованием технологии CUDA;
•	навыками разработки, компиляции и отладки параллельных программ с использованием технологии OpenCL;
•	навыками разработки, компиляции, профилирования и отладки параллельных программ с использованием технологии MPL
При подготовке этого издания были существенно переработаны темы 1—4 и добавлена тема 5, посвященная программированию гетерогенных компьютеров.
Тема!
АРХИТЕКТУРА ПАРАЛЛЕЛЬНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ
Параллельные вычислительные системы (ПВС) являются стратегическим направлением развития средств вычислительной техники [1]. Это обусловлено следующими причинами.
С одной стороны, всегда существовали, существуют сейчас, и, по всей вероятности, будут существовать всегда такие задачи, которые невозможно решить за приемлемое время на одиночном вычислительном устройстве — процессоре или ядре многоядерного процессора [2]. Необходимость решать большие задачи постоянно стимулировала вложение средств и усилий в создание все более быстродействующих процессоров. При этом рост производительности обеспечивался как за счет технологического совершенствования, так и путем реализации разнообразных видов внутреннего или скрытого параллелизма: конвейеризации, кэширования памяти, суперскалярности, векторизации [3]. Тем не менее, с точки зрения программиста, использующего обычные широко распространенные языки программирования, такие процессоры все равно являются устройствами, строго последовательно реализующими заданные шаги алгоритма решения задачи.
С другой стороны, вычислительные возможности одиночного последовательного устройства оказались существенно ограничены физическими законами. Невозможно беспредельно увеличивать тактовую частоту, плотность размещения транзисторов на кристалле и другие параметры, от которых зависит производительность. Один процессор, а точнее — одно ядро процессора, к настоящему моменту времени (первая четверть XXI в.) способен выполнить примерно 109—1010 операций с плавающей точкой в секунду (от одного до 10 гигафлопс). Ожидать резкого увеличения производительности одиночного вычислителя не приходится. Однако такой производительно
7
сти заведомо недостаточно для решения ряда больших задач в таких областях, как прогнозирование погоды, анализ генома, ядерная энергетика, моделирование самолета/ракеты, сейсмогеофизика и т. д.
В силу этого наряду с совершенствованием одиночных последовательных вычислителей практически постоянно исследовались и создавались коллективы таких вычислителей — параллельные вычислительные системы. Общий объем вычислений, необходимых для решения задачи, распределяется между процессорами, входящими в коллектив, таким образом, чтобы каждый из них мог выполнять свою часть работы одновременно с остальными. Параллельность вычислений позволяет получить эффект ускорения, т. е. уменьшения затрат времени на решение задачи по сравнению с одиночным последовательным вычислителем. Термин «коллектив» подразумевает, что его составные части могут взаимодействовать друг с другом в процессе решения единой задачи. Программа решения задачи для коллектива вычислителей уже не может быть полностью выражена на каком-либо строго последовательном языке программирования, в ней должен быть определен именно явный параллелизм работы составных частей системы и описано их взаимодействие.
Архитектура параллельных систем с самого начала их создания и применения развивалась в самых различных направлениях, исследовалось большое количество вариантов реализации параллелизма. Огромное разнообразие архитектур вычислительных систем породило естественное желание ввести для них какую-то классификацию. Эта классификация должна позволять однозначно относить ту или иную вычислительную систему к некоторому классу, который, в свою очередь, должен достаточно полно характеризовать ее основные свойства. Попыток классификации предпринималось множество. Наиболее популярная классификация была предложена М. Флинном в конце 1960-х гг. [3]. Она использует понятия двух потоков: команд и данных. На основе соотношения количеств этих потоков выделяется четыре основных класса архитектур.
1.	ОКОД (Single Instruction Single Data, SISD) — единственный поток команд и единственный поток данных. Этим определяется классическая машина фон Неймана. К классу ОКОД относятся все однопроцессорные (одноядерные) компьютеры.
8
2.	ОКМД (Single Instruction Multiple Data, SIMD) — единственный поток команд и множественный поток данных. Типичными представителями этого класса являются матричные компьютеры, в которых все процессорные элементы выполняют в едином темпе одну и ту же операцию, задаваемую единственным устройством управления. Каждый процессор обрабатывает свои собственные данные. В настоящее время к этому классу систем относятся графические процессоры, появившиеся в результате развития видеоадаптеров персональных компьютеров путем переноса в них функций обработки пикселов изображения из центрального процессора.
3.	МКОД (Multiple Instruction Single Date, MISD) — множественный поток команд и единственный поток данных. Представителями такой архитектуры можно считать высоконадежные вычислительные системы с дублированием процессов обработки одних и тех же данных на нескольких вычислительных устройствах и мажоритарным выбором результатов решения задачи.
4.	МКМД (Multiple Instruction Multiple Date, MIMD) — множественный поток команд и множественный поток данных. Входящие в систему вычислительные устройства каждое в своем темпе обрабатывают собственные данные в соответствии с собственной (но, как правило, одинаковой для всех) программой. Это чрезвычайно обширный класс параллельных вычислительных систем, который принято делить на два подкласса: мультипроцессоры и мультикомпьютеры.
Мультипроцессоры или системы с общей памятью состоят из нескольких (как правило, одинаковых) процессоров, имеющих равноприоритетный доступ к общей памяти с единым адресным пространством (ОП-МКМД или SM-MIMD). К этому подклассу относятся [3,4]:
•	векторно-конвейерные компьютеры;
•	SMP-системы (symmetric multi processors');
•	UMA-системы (uniform memory access);
•	NUMA-системы (non-uniform memory access);
•	ccNUMA-системы (cache coherent NUMA);
•	многоядерные процессоры.
К сожалению, производительность систем этого подкласса нелинейно изменяется при увеличении количества процессоров и практически перестает расти, начиная с 16—20 (исключение составляют ccNUMA-системы, для которых максималь
9
ное количество процессоров может быть порядка нескольких сотен) [3].
В последние годы стремительно развивающимся направлением развития ПВС стал подкласс мультикомпьютеров — систем с распределенной памятью (РП-МКМД или DM-MIMD). Идея построения вычислительных систем данного подкласса очень проста: некоторое количество самостоятельных вычислительных узлов объединяются друг с другом высокоскоростной коммуникационной средой. Каждый вычислительный узел имеет один или несколько процессоров и свою собственную локальную память, разделяемую этими процессорами. Распределенность памяти означает то, что каждый процессор имеет непосредственный доступ только к локальной памяти своего узла. Доступ к данным, расположенным в памяти других узлов, выполняется значительно дольше и существенно более сложными способами.
Компоненты параллельной системы с распределенной памятью могут проектироваться специально для одного образца или небольшой серии. Такие системы называют массивно-параллельными процессорами (МПП). Основным преимуществом таких систем является значительно большая масштабируемость по сравнению с системами с общей памятью. В зависимости от класса решаемых задач и бюджета может быть реализована система с числом узлов от нескольких десятков до нескольких сотен тысяч. Однако сверхвысокая стоимость массивно-параллельных процессоров, обусловленная уникальностью их компонентов, не позволяет применять их во многих областях, нуждающихся в системах высокой производительности. Это привело к развитию вычислительных кластеров — аналогов МПП, но собранных не из уникальных компонентов, а из тех, которые можно свободно приобрести на рынке.
Технологической основой развития кластеризации стали широкодоступные и относительно недорогие микропроцессоры и коммуникационные (сетевые) технологии, появившиеся в 1990-х гг. Вычислительный кластер представляет собой совокупность вычислительных узлов (от нескольких десятков до сотен тысяч), управляющей системы и системы хранения данных. Структура (более-менее типичная) вычислительного кластера высокой производительности или массивно-параллельного процессора показана на рис. 1.1. Вычислительные узлы обычно объединяются как минимум двумя, обычно независимыми,
ю
сетями передачи данных. Управляющая сеть используется для распределения программ и исходных данных, запуска параллельных программ на выполнение и сбора результатов вычислений. Основная коммуникационная сеть (как правило, значительно более производительная по сравнению с управляющей) обеспечивает обмен данными между программами, исполняемыми на узлах в процессе вычислений. В самых высокопроизводительных системах к ним обычно добавляются сеть синхронизации и сеть тестирования. Управляющая система имеет выход в Интернет для доступа удаленных пользователей к ресурсам кластера.
В настоящее время в списке пятисот самых высокопроизводительных систем мира свыше 450 — это кластеры. Все остальные относятся к подклассу массивно-параллельных процессоров. Дополнительные сведения об архитектуре параллельных вычислительных систем можно найти в работах [4—10].
Рис. 1.1. Вычислительный кластер или массивно-параллельный процессор
Тема 2
ВВЕДЕНИЕ В ПАРАЛЛЕЛЬНОЕ ПРОГРАММИРОВАНИЕ
Если некоторую задачу не удается решить за приемлемое время на одиночном последовательном вычислителе, то приходится разрабатывать программу ее решения на параллельной вычислительной системе. В некоторых случаях большую задачу удается разбить на несколько более мелких независимых друг от друга подзадач. Тогда для решения каждой из них можно по-прежнему писать последовательные программы на обычных языках программирования и получать требуемый эффект ускорения просто за счет одновременного запуска этих программ на разных, возможно даже физически несвязанных вычислительных устройствах. Это называется параллелизмом задач или внешним параллелизмом. В связи с тем, что каждая программа, очевидно, обрабатывает свои собственные данные, такой способ обозначается как МПМД — много программ, много потоков данных (в английской литературе — MPMD) [3].
Однако реально существуют ситуации, когда разбиение задачи на несвязанные подзадачи либо невозможно, либо не дает желаемого ускорения времени ее решения. В этом случае для более быстрого решения всей задачи в целом или ее отдельных подзадач на параллельной вычислительной системе приходится использовать так называемый внутренний параллелизм алгоритма. Под внутренним параллелизмом алгоритма понимается потенциальная возможность одновременного выполнения некоторых составляющих его операций.
С теоретической точки зрения внутренний параллелизм алгоритма или реализующей его программы может быть выражен графом информационных зависимостей между операциями, а точнее — его ярусно-параллельной формой [1,2]. В этой форме все операции, зависящие только от значений исходных данных, располагаются в виде вершин на верхнем
12
ярусе. Те операции, которые используют результаты, вырабатываемые на верхнем ярусе, размещаются на следующем ярусе графа, причем зависящие друг от друга вершины соединяются дугами. Этот процесс продолжается до тех пор, пока в граф не будут включены все операции алгоритма/программы.
Преобразование последовательного алгоритма/программы в ярусно-параллельную форму и последующее исследование свойств и характеристик этого графа могут быть основой для составления расписания одновременного выполнения всех не зависящих друг от друга операций. Такой подход называется мелкозернистым распараллеливанием, теоретически он позволяет получить максимальное ускорение. Однако из-за большой трудоемкости к мелкозернистому распараллеливанию обычно прибегают в том случае, когда менее формализованные способы не позволяют добиться нужного результата.
Существует другой подход, позволяющий на практике и при ручном применении получить вполне приемлемые значения ускорения процесса вычислений за счет так называемого крупнозернистого [5] (или крупноблочного [6]) распараллеливания. Известно, что более 90 % (как правило, значительно более, т. е. 99.99(9) %) из времени исполнения программы приходится на 10 % (обычно — существенно меньше) ее кода: так называемое соотношение 90:10. Такими участками кода являются циклы и гнезда вложенных друг в друга циклов. Именно этими участками текста программы и следует заниматься при ее крупнозернистом распараллеливании.
Нециклические (чисто последовательные) участки кода программы распараллеливать никакого смысла не имеет. Причина состоит в том, что при современной производительности одиночного процессора/ядра затраты времени на выполнение нециклической последовательности операций длиной, например, в миллион машинных команд, составят меньше одной миллисекунды. Представить себе существование такой или более длинной реальной нециклической программы очень сложно — их просто нет. Практически любая программа содержит циклы, а те, которые нуждаются в распараллеливании, т. е. имеют большое время выполнения, безусловно, являются циклическими. Поэтому предметом дальнейшего рассмотрения при распараллеливании являются циклы и гнезда вложенных друг в друга циклов.
13
2.1.	Распараллеливание циклов
Теперь в качестве простой операции имеет смысл рассматривать тело некоторого цикла полностью, каким бы сложным оно ни было. Если между какими-либо двумя итерациями цикла существует информационная зависимость (значения, вырабатываемые при выполнении одной итерации, используются в другой), то выполнять эти итерации одновременно невозможно.
Приведем простой пример такого цикла. Пусть необходимо вычислить сумму элементов массива агг. Типичный фрагмент последовательной программы на С-подобном языке может выглядеть так:
double arr[ size ];
// формирование значений элемента массива double sum = 0;
for (int i = 0; i < size; i++) sum += arr[ i ];
Каждая итерация данного цикла, кроме самой первой, использует значение переменной sum, вычисленное в результате выполнения всех предыдущих итераций, что и называется информационной зависимостью итераций цикла. Этот участок программы впрямую не может быть распараллелен, однако существуют приемы и способы, позволяющие преобразовать некоторые циклы с информационной зависимостью к виду, допускающему параллельное выполнение итераций. Такие приемы могут существенно зависеть от предполагаемого способа исполнения параллельной программы, т. е. от технологии подготовки и runtime-поддержки, которая, в свою очередь определяется архитектурой и классом/подклассом параллельной вычислительной системы. Некоторые из них, реализованные в соответствующих технологиях поддержки параллельных вычислений, будут рассмотрены в последующих темах пособия.
Таким образом, выполняться параллельно могут только такие циклы, итерации которых информационно не зависят друг от друга. Код итерации такого цикла и является «крупным зерном», т. е. единичной операцией, целиком назначаемой на выполнение одному вычислительному устройству. При этом, как правило, чем больший объем вычислений выполняется одним таким «зерном», тем более эффективным получается распа
14
раллеливание. В результате применения крупнозернистого способа распараллеливания последовательная программа преобразуется в такую форму, которая будет выполняться одновременно несколькими вычислительными узлами параллельной системы. Каждый узел, выполняя одну и ту же программу, обрабатывает свои собственные данные. Это называется режимом ОПМД — одна программа, много потоков данных (SPMD в английском варианте).
До настоящего момента по существу предполагалось, что целевая вычислительная система, которая будет исполнять распараллеленную программу, имеет настолько много вычислительных узлов, что каждому узлу достанется ровно по одному зерну (итерации цикла) на исполнение. Другими словами, не рассматривался вопрос о распределении множества итераций цикла между физическими вычислителями. Реальная ситуация может быть совершенно другой. Количество узлов системы может быть меньше или даже много меньше количества итераций цикла. Разработка параллельной программы не может быть выполнена без решения вопроса о способе распределения итераций каждого распараллеливаемого цикла между процессорами системы. Пусть в программе есть цикл, между итерациями которого нет информационной зависимости, например:
for(i =0; i < N; i += 1)
arr[i] = <выражение, не содержащее arr[k] (k! = i)>;
В дальнейшем для краткости вместо правой части оператора присваивания <выражение, не содержащее arr[k] (k! = i)> будем использовать обозначение/(z,...).
Итерации этого цикла при исполнении программы на одном вычислителе выполняются строго последовательно для значений индекса цикла 0, 1, 2, ..., N - 1.
При исполнении программы на параллельной системе, содержащей и вычислительных узлов (процессоров) эти N итераций могут быть по-разному распределены между ветвями. Ветвью будем называть экземпляр параллельной программы, исполняемый одним узлом. Заметим, что в принципе возможна ситуация, когда один узел выполняет не одну ветвь, а несколько, в пределе — все ветви (например — в процессе отладки).
Самыми часто используемыми являются следующие способы распределения: блочное, циклическое, блочно-цикличе
15
ское, динамическое с фиксированным и динамическое с переменным уменьшающимся размером блока итераций.
Блочное распределение осуществляется по N/n итераций на узел по порядку возрастания индексов цикла (если N не делится нацело на п, то разным узлам достанется разное количество итераций). При этом вычислительные узлы фактически будут выполнять следующие фрагменты программы (в этих фрагментах через NO, N1,..., Nk обозначены границы значений индекса цикла для соответствующих узлов):
for(i = 0; i < N0; i += 1) // доля узла 0, N0 = N / п arr[i] = f(i, ...);
for(i = N0; i<Nl; i += 1) // доля узла 1, N1 = 2*N I n arr[i] = f(i, ...);
for(i = Nk; i < N; i += 1) // Nk = (n-l)*N/n arr[i] = f(i, ...);
На рис. 2.1 в графическом представлении изображено блочное распределение пространства итераций между п узлами.
NO-2	NO —1
N1 - 2	N1 - 1
N-1
Рис 2.1. Блочное распределение
Этот способ очень просто реализуется, не требует дополнительных накладных расходов на управление вычислениями и довольно равномерно распределяет вычислительную нагрузку между узлами при условии одинаковости затрат времени на выполнение каждой итерации. Чем равномернее распределен весь объем вычислительной работы между параллельно функционирующими узлами, тем выше коэффициент ускорения решения задачи. Однако равномерность загрузки может быть сильно нарушена, если длительности выполнения блоков итераций существенно различаются (в том числе потому, например, что узлы системы имеют разную производительность или на некоторых из них выполняются прочие конкурирующие процессы). Это справедливо и для
16
остальных статических способов (циклического и блочноциклического), поэтому существуют динамические способы. Применимость того или иного способа распределения зависит от многих факторов (в том числе от архитектуры вычислительной системы).
При циклическом распределении (тоже по N/n итераций) фрагменты, исполняемые разными узлами, выглядят так:
for(i = 0; i<N; i +=n) 11 ppim узла 0 arr[i] = f(i, ...);
for(i = 1; i<N; i += n) 11 ppnn узла 1
arr[i] = f(i, ...);
for(i = n - 2; i < N; i += n) // доля узла n - 1 arr[i] = f(i, ...);
Графическое представление циклического распределения изображено на рис. 2.2.						
Узел 0:	0	п	2п	• ••	(т - 1)п	тп
	...					
Узел к:	к	п + к	2п + к	...	(т - 1)п + к	N-1
	...					
Узел п -1:	п-1	2п-1	Зп-1	...	тп -1	
Puc 2.2. Циклическое распределение
Блочно-циклическое распределение (это циклическое распределение, но каждый узел исполняет несколько блоков итераций):
for(k =0; k < N; к += р)
for(i = к * n; i < (к + 1) * n; i += 1) // доля узла 0 arr[i] = f(i,..);
for(k =0; к < N; к += р)
for(i = (к + 1) * n; i < (к + 2) * n; i += 1) // доля
узла 1
arr[i] = f(i,..);
Графическое представление блочно-циклического распределения изображено на рис. 2.3.
17
Узел 0:	0	1	2	(и - 1)р + 3
Узел 1:	3	4	5	(и - 1)р + 6
	...			
Узел п -1:	(и - 1)р	(п-1)р + 1	(и - 1)р + 2	2(п-1)р
Рис 2.3. Блочно-циклическое распределение
N-1
Как уже отмечалось, для всех статических методов характерна возможность неодинаковых затрат времени параллельно работающих узлов на выполнение заранее заданного фиксированного объема работы. Следствием этого является возможность уменьшения коэффициента ускорения по сравнению с потенциально достижимым.
Динамические методы распределения ориентированы на сглаживание этой неравномерности за счет того, что блоки итераций размером (N / и) / к (при к > 1) выделяются узлам по мере их готовности. Использование этих методов предполагает наличие некоего единственного «центра», управляющего списком невыполненных блоков итераций. Каждый узел в самом начале выполнения распараллеленного цикла или после выполнения ранее полученного блока обращается к этому центру, запрашивая следующий блок. Взаимодействия узлов с центром распределения блоков являются, по существу, накладными расходами, суммарный объем которых должен быть значительно меньше затрат времени на выполнение блоков итераций. Чем меньше размер блока итераций, тем меньше потенциальная неравномерность загрузки, но тем больше накладные расходы. Поэтому при динамическом распределении блоками итераций фиксированного размера должен быть найден оптимальный размер блока, при котором малы накладные расходы, и не проявляется неравномерность загрузки узлов. Иногда наилучшим решением будет уменьшение размера блока итераций вплоть до единицы в процессе вычислений, т. е. реализация динамического распределения с переменным размером блока. При этом можно ожидать сравнительно небольшой доли накладных расходов вследствие большого размера блоков итераций в начале выполнения цикла и сглаживания неравномерности загрузки узлов за счет именно динамического распределения, а также малого размера блоков в конце вычислений. Графическое пред
18
ставление, которым иллюстрировались статические методы распределения, для динамических методов сформировать в общем виде затруднительно, поэтому здесь оно не приводится.
Алгоритмы единого центра динамического распределения блоков итераций отнюдь не просты и не могут быть рекомендованы для самостоятельной реализации при распараллеливании каждого цикла программы на системах, не имеющих общей памяти. Для систем с общей памятью и статические, и динамические способы распределения реализованы в технологии ОрепМР, рассматриваемой в теме 3.
2.2.	Распараллеливание гнезд циклов
Для больших задач типичной является ситуация, когда основной объем вычислений определяется совокупностью вложенных друг в друга циклов (гнездом циклов):
for(i = 0; i < N; i += 1)
for(j = 0; j < M; j += 1)
for(k = R - 1; k >= 0 ; k -= 1)
<тело внутреннего цикла>
Распределение итераций только одного внешнего цикла может оказаться невозможным из-за их взаимной информационной зависимости или нецелесообразным вследствие малого их количества и, как следствие, неравномерности загрузки и невысокого коэффициента ускорения. Поэтому может оказаться необходимым распределение всего пространства итераций гнезда вложенных циклов. Таким пространством называют множество целочисленных векторов {[О, О, R - 1,...], [О, О, R -2], ...}, координаты которых задаются значениями индексных переменных всех циклов данного гнезда. Задача распараллеливания в этом случае сводится к разбиению этого множества векторов на подмножества таким образом, чтобы итерации этих подмножеств могли быть выполнены одновременно.
Существуют следующие методы разбиения: координат, гиперплоскостей, параллелепипедов, пирамид и более сложные.
Метод координат заключается в том, что пространство итераций фрагмента разбивается на гиперплоскости, перпендикулярные одной из координатных осей.
19
Пусть, например, есть гнездо циклов:
for(i = 1; i < N; i += 1)
for(j = 0; j < M; j += 1)
a[i][j] =a[i - l][j] + a[i][j];
Разбиение пространства итераций по измерению i (гиперплоскостями, перпендикулярными оси i) приводит к разрыву информационных зависимостей. Следовательно, такое распараллеливание неприемлемо. Однако возможно применение метода координат с разбиением пространства итераций гиперплоскостями, перпендикулярными оси j. Все такие операции назначаются для выполнения на одном узле системы.
При использовании метода гиперплоскостей пространство итераций размерности п разбивается на подпространства размерности от 1 до и - 1 таким образом, чтобы все операции, соответствующие точкам одного подпространства, могли выполняться параллельно. Каждое такое подпространство есть гиперплоскость, не обязательно перпендикулярная какой-либо оси координат.
Например, для такого гнезда:
for(i =1; i < N; i += 1)
for(j = 1; j < M; j += 1)
a[i][j] =a[i - l][j] + a[i][j - 1];
метод координат прямо не может быть применен. Однако существуют гиперплоскости, удовлетворяющие условию i + j = = const, и проходящие через вершины, не связанные информационными зависимостями. При этом параллельно работающим узлам следует назначить вычисления тела цикла для каждой из гиперплоскостей (совокупностей итераций, для которых i+j = const).
Метод параллелепипедов является логическим развитием двух предыдущих методов и заключается в разбиении пространства итераций на n-мерные параллелепипеды с использованием ортогональных гиперплоскостей, весь объем каждого из которых назначается отдельному процессору.
Существуют и другие, более сложные методы распределения пространства итераций, так же основанные на том, что параллельно выполняемые итерации не должны информационно зависеть друг от друга [7—10].
20
2.3.	Преобразования гнезд циклов
Возможны ситуации, когда гнездо циклов исходной последовательной программы не может быть распараллелено никаким из описанных методов. Тогда можно попытаться применить одно или несколько эквивалентных преобразований гнезда циклов для изменения характера информационных зависимостей между итерациями. К таким преобразованиям относятся следующие.
Перестановка циклов: внешний и внутренний циклы меняются местами (если это возможно). Например, фрагмент:
for(i = 1; i < N; i += 1)
for(j = 1; j < M; j += 1)
a[i][j] = a[i - l][j] + a[i][j - 1];
заменяется эквивалентным фрагментом:
for(j = 1; j < M; j += 1)
for(i = 1; i < N; i += 1)
a[i][j] =a [i - l][j] + a[i][j - 1];
Распределение цикла: один цикл расщепляется на два цикла или более. В примере из перестановки циклов то же самое гнездо переписано так, что один внутренний цикл заменен на два. Итерации внутреннего цикла теперь будут выполняться в другой последовательности, но общий объем вычислительной работы остался неизменным:
for(i = 1; i < N; i += 1) {
for(j = 1; j < M ; j += 2)
a[i][j] = a[i - l][j] + a[i][j - 1];
for(j = 2; j < M ; j += 2)
a[i][j] = a[i - l][j] + a[i][j - 1];
}
Слияние циклов: два цикла объединяются в один. В том же примере исходное гнездо двух циклов превращается в один цикл (общий объем вычислительной работы остался таким же):
for(i = 1, j = 1; (i < N) && (j < M);
i += ++j/M, j = (j >= M? 1 : j))
a[i][j] = a[i - l][j] + a[i][j - 1];
21
Развертка цикла: тело цикла записывается два раза (или большее количество) при соответствующих изменениях способов модификации индексных переменных. Общий объем вычислений остается тем же самым, но количество итераций уменьшается:
for(i = 1; i < N; i + = 1)
for(j = 1; j < M; j += 2){
a[i][j] = a[i - l][j] + a[i][j - 1];
a[i][j + 1] = a[i - 1][j +1] + a[i][j];
}
Расширение скаляра: одиночная переменная, значение которой вычисляется предыдущими итерациями и используется в последующих итерациях, заменяется элементом массива. Для примера рассмотрим фрагмент программы умножения матриц:
for(i = 0; i < N; i += 1)
for(j = 0; j < N; j += 1){
sum = 0.0;
for(k =0; k < N; k += 1) sum += a[i][k] *b[k][j]; c[i][j] = sum;
}
Чтобы ликвидировать информационные зависимости итераций самого внутреннего цикла друг от друга, используем целевой элемент массива вместо промежуточной переменной:
for(i =0; i < N; i += 1)
for(j = 0; j < N; j += 1){
c[i][j] =0.0;
for(k =0; k < N; k += 1) c[i][j] += a[i][k] *b[k][j];
}
В результате этого преобразования исчезла зависимость итераций всех циклов гнезда друг от друга и соответственно появилась возможность распараллелить их.
Разделение на полосы: один цикл разбивается на два вложенных цикла. В предыдущем примере самый внешний цикл разобьем на два, добавив новую индексную переменную ind и шаг разбиения step (если N нацело делится на step, то все полосы будут иметь равную ширину, иначе последняя будет меньше):
22
for(ind = 0; ind<N;ind += step)
for(i = ind; i < ind + step; i += 1)
for(j = 0; j < N; j += 1){
c[i][j] =0.0;
for(k =0; к < N; к += 1) c[i][j] += a[i][k] *b[k][j];
}
Обычно разделение на полосы используется как часть разделения на блоки (см. следующий способ), развертки и сжатия.
Разделение на блоки: область итераций разбивается на прямоугольные блоки. Целью такого преобразования может быть согласование с размером области основной памяти, перемещаемой в кэш процессора при обращении к любому ее элементу.
Рассмотрим гнездо циклов транспонирования матрицы:
for(i =0; i < N; i += 1)
for(j = 0; j < N; j += 1)
b[i][j] = a [j][i];
Пусть известна ширина линии кэша, назовем ее width (это значит, что при обращении к элементу a [j] [z] в кэш переносится width элементов строки матрицы а, которые в ближайших итерациях не будут использоваться и будут вытеснены). Преобразуем гнездо из двух циклов в гнездо, содержащее четыре цикла:
for(i =0; i < N; i += width)
for(j = 0; j < N; j += width)
for(ii =0; ii < min(N, i + width - 1); ii += 1)
for(jj = 0; jj< min(N, j + width - 1); jj += 1) b[ii][jj] = a [jj][ii];
Теперь два внутренних цикла работают с фрагментом матрицы размером (width х width). Элементы массива а, перенесенные в кэш при обращении к самому левому элементу фрагмента, будут выбираться при исполнении ближайших итераций внутреннего цикла без обращений к основной памяти. Поэтому транспонирование матрицы будет выполняться значительно быстрее.
Развертка и сжатие: комбинируются перестановка циклов, разделение на полосы и развертка.
23
Перекос цикла: выполняется модификация границ изменения индексов циклов таким образом, чтобы сделать явной параллельность фронта волны изменения элементов массивов. Например, при решении дифференциальных уравнений в частных производных обычно используется цикл
for(i = 0; i < N; i += 1)
for(j = 0; j < N; j += 1)
a[i][j] = (a[i - l][j] + a[i][j - 1] + a[i + 1] [j]
+ a[i] [j +1]) / 4;
После внесения такой модификации (перекоса цикла):
for(i =0; i < N; i += 1)
for(j = i; j < i + N; j += 1)
a[i][j - i] = (a[i - 1][j - i] + a[i][j - i - 1]
+ a[i +1] [j - i] + a[i][j - i + 1]) / 4;
явно обнаружится возможность распараллеливания по столбцам, для чего надо будет просто переставить циклы по i и по j.
2.4.	Проблемы разработки параллельных программ
Существуют различные подходы к разработке параллельных программ, ориентированные на разные архитектуры высокопроизводительных вычислительных систем и различные инструментальные средства.
1.	Модель общей памяти. Согласно этой модели параллельно исполняемые ветви одной программы (потоки управления) обращаются к общей памяти для выполнения операций считы-вания/записи. В рамках этой модели не требуется описывать обмен данными между задачами в явном виде. Это упрощает программирование. Вместе с тем особое внимание приходится уделять соблюдению детерминизма и таким явлениям, как «гонки за данными», блокировки и т. д.
2.	Модель передачи сообщений, характеризуемая тем, что программа порождает несколько задач. Каждой задаче присваивается свой уникальный идентификатор. Взаимодействие задач осуществляется посредством отправки и приема сообщений. Новые задачи могут создаваться во время выполнения
24
параллельной программы, несколько задач могут выполняться на одном процессоре.
3.	Модель параллелизма данных, для которой характерно следующее:
•	одна и та же операция применяется к множеству элементов структуры данных;
•	гранулярность вычислений невелика;
•	программист должен указать транслятору, как данные следует распределить между параллельными ветвями.
Независимо от используемой модели основная сложность при проектировании параллельных программ состоит в том, чтобы обеспечить правильную последовательность взаимодействий между различными вычислительными процессами, а также координацию использования ресурсов, разделяемых между процессами. При разработке параллельных программ обычно возникает целый ряд специфических по сравнению с последовательным программированием проблем:
•	необходимость управления выполнением множества параллельно протекающих вычислительных процессов;
•	необходимость заботиться о масштабируемости программы (масштабируемость — это свойство вычислительной системы, компьютерной сети или совокупности процессов, которое заключается в сохранении способности эффективно работать при возрастании количества узлов);
•	необходимость заботиться о сбалансированной загрузке всех используемых вычислительных узлов;
•	необходимость обеспечения межпроцессных взаимодействий, контроль их корректности, эффективности и степени влияния на время вычислений;
•	возможность возникновения тупиковых ситуаций или взаимных блокировок при управлении доступом к общим ресурсам и, как следствие, необходимость их предотвращения;
•	нелокальный и динамический характер ошибок, существенным образом усложняющий отладку программы;
•	утрата детерминизма (предсказуемого поведения программы в процессе ее исполнения). Параллельная программа может давать разные результаты при повторных запусках даже без модификации кода и исходных данных. Это является следствием так называемых «гонок за данными» — одновременного и асинхронного доступа разных процессов или потоков к общим переменным.
25
Наиболее специфическими проблемами, возникающими при разработке параллельных программ в модели с общей памятью, является возможность возникновения:
•	гонок данных (data race);
•	блокировок (deadlock).
Гонки данных. Они являются следствием таких информационных зависимостей, когда несколько вычислительных процессов нескоординировано модифицируют содержимое одной и той же области памяти (переменной). Наличие возможности возникновения гонок данных в программе далеко не всегда является очевидным. Пусть, например, два (или больше) процесса время от времени выполняют, казалось бы, очень простую операцию увеличения на единицу значения общей переменной:
count += 1;
Эта операция, представляющаяся неделимой, на самом деле состоит как минимум из трех действий.
1.	Чтение из памяти в регистр процессора/ядра значения переменной count.
2.	Увеличение содержимого регистра на 1.
3.	Запись значения из регистра в память.
Если в промежутке времени между операциями чтения и записи, выполняемыми одним процессом, какой-либо другой процесс прочитает значение переменной в свой регистр, то после того как оба процесса выполнят запись, переменная count будет иметь значение на единицу большее, чем начальное. Но это неверно, поскольку два процесса добавляли по единице, и значение count должно было увеличиться на 2. Если же во время работы параллельной программы в промежутке между чтениями и записями, выполняемыми любым процессом, никакой другой процесс не будет обращаться к переменной count, то ее значение будет формироваться правильно.
Первая ситуация и называется «гонками данных». Параллельная программа может выдавать совсем не те результаты, которые выдала бы последовательная программа. Такая параллельная программа может вести себя по-разному при разных ее запусках в зависимости от временных соотношений моментов доступа ее процессов к общим данным.
26
Имеются два способа борьбы с гонками данных:
•	использование только локальных для каждого процесса, а не разделяемых переменных. Однако полностью без разделяемых переменных в параллельных программах обойтись невозможно, если их процессы действительно выполняют вычисления, направленные на решение одной задачи;
•	управление доступом к разделяемым переменным с помощью различных средств синхронизации таким образом, чтобы цепочка действий «чтение — модификация — запись» одного процесса не могла быть прервана другим процессом. Средства синхронизации могут быть реализованы с помощью атомарных операций, критических секций, семафоров, мьютексов и других механизмов. Неосторожное или неграмотное применение некоторых из этих средств может привести к возникновению второй серьезной проблемы — блокировкам.
Блокировка. Она возникает, если процесс ожидает выполнение условия, которое не может быть выполнено. Например, пусть с помощью какого-либо из средств синхронизации обеспечивается корректный доступ к модификации двух общих переменных varl и var2. Пусть теперь один процесс получил монопольный доступ к переменной varl и пытается получить такой же доступ к переменной var2, но ее уже захватил другой процесс, который, в свою очередь, пытается получить в свое пользование переменную varl. Если синхронизация доступа к общей переменной сопровождается приостановкой процесса до момента, когда доступ может быть предоставлен, то оба процесса будут остановлены (взаимно заблокированы) и не смогут предпринять никаких действий по преодолению сложившейся ситуации. Разрешить ее возможно только извне путем остановки параллельной программы, т. е. принудительного прекращения решения задачи.
Активной блокировкой называют ситуацию, когда процессы не приостанавливаются (занимают вычислительный ресурс), но и не производят полезных вычислений, непрерывно пытаясь получить доступ к захваченным другими процессами ресурсам. Это возможно в случае использования некоторых средств синхронизации, таких как omp_test_lock() в технологии ОрепМР. В этом случае в программе можно предусмотреть возникновение блокировок и их разрешение, но это требует значительных усилий. Особенно тогда, когда общих ресурсов не два, а много, и так же много нуждающихся в них процессов.
27
Для предотвращения блокировок можно рекомендовать, например:
•	пронумеровать все общие ресурсы;
•	захватывать их строго в порядке возрастания номеров.
Тогда вышеописанной ситуации блокировки возникнуть не может. К сожалению, эта технология не может быть с успехом применена во всех приложениях. Например, в условиях высоконагруженных событийно-управляемых вычислений последовательность возникновения потребности в общих ресурсах обычно совершенно непредсказуема, а необходимость захватывать ненужные ресурсы просто потому, что такова дисциплина, приводит к падению производительности.
Программирование для параллельной вычислительной системы самым существенным образом зависит от ее архитектуры и от используемой технологии поддержки параллельных вычислений. Под технологией поддержки здесь понимается совокупность программных средств, используемых при подготовке и для исполнения параллельной программы. Для разных классов/подклассов параллельных систем существуют ориентированные на их специфику специально разработанные эффективные программные средства поддержки разработки и исполнения пользовательских программ.
В этом пособии освещаются следующие технологии:
•	ОрепМР — для систем класса МКМД/ОПМД с общей памятью (многоядерные процессоры);
•	CUDA — для систем класса ОКМД (графические процессоры компании NVIDIA);
•	OpenCL—для гетерогенных систем класса МКМД + ОКМД (многоядерные процессоры с графическими процессорами или вычислительными акселераторами);
•	MPI — для систем класса МКМД/ОПМД с распределенной памятью (вычислительные кластеры, массивно-параллельные процессоры).
ТемаЗ OPENMP: ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ НА СИСТЕМАХ С ОБЩЕЙ ПАМЯТЬЮ
Технология ОрепМР [11,12] в настоящее время является одним из наиболее популярных средств параллельного программирования для компьютеров с общей памятью. Она базируется на традиционных последовательных языках программирования и использовании комментариев специального вида в качестве указаний компилятору о том, как следует преобразовать оператор (блок операторов) при создании параллельной программы. За основу берется последовательная программа, а для создания ее параллельной версии пользователю предоставляется набор директив, функций и переменных окружения. Технология ОрепМР нацелена на то, чтобы пользователь имел один и тот же вариант текста программы как для параллельного, так и для последовательного режима выполнения.
Важным достоинством технологии ОрепМР является возможность реализации так называемого инкрементального программирования. При этом подходе программист постепенно находит в программе участки, содержащие ресурс параллелизма, с помощью предоставляемых ею механизмов делает их параллельными, а затем переходит к анализу следующих участков. Таким образом, распараллеленные участки постепенно охватывают всё большую часть программы. Этот подход значительно облегчает процесс адаптации последовательных программ к параллельным компьютерам, а также отладку и оптимизацию программ.
3.1.	Модель параллельной программы ОрепМР
Модель исполнения параллельной программы, подготовленной с помощью технологии ОрепМР, можно сформулировать следующим образом [13]:
29
•	при запуске программы создается главный поток, выполняющий впоследствии все последовательные области программы;
•	программа содержит набор последовательных и параллельных областей (или секций или регионов);
•	при входе в параллельную область главным потоком выполняется операция Fork, порождающая совокупность подчиненных потоков. Каждый поток имеет свой уникальный числовой идентификатор (главному потоку соответствует 0);
•	потоки одного региона выполняют один и тот программный код (с точностью до возможных разветвлений), но каждый в идеале обрабатывает свои собственные данные. Если программный код региона представляет собой распараллеливаемый цикл, то разные потоки выполняют разные его итерации. Это распределение вычислительной работы между потоками и позволяет добиваться ускорения выполнения программы по сравнению с ее последовательным вариантом;
•	при выходе из параллельной области всеми потоками выполняется операция Join. Завершается выполнение всех пото
ков, кроме главного.
Создание и исполнение двух параллельных областей иллюстрирует рис. 3.1. В первой области все потоки приходят к моменту выхода из нее практически одновременно, во второй — в разные моменты времени (как правило, это более соответствует реальной динамике выполнения потоков).
Основной (master) поток:
Параллельные регионы
Рис 3.1. Создание двух параллельных регионов
Модель памяти, используемая в ОрепМР. Все переменные, объявленные в программе вне фрагментов, выполняемых в качестве параллельных регионов, по умолчанию относятся к классу разделяемых (shared), т. е. общедоступных для всех потоков.
30
Если программа вошла в параллельный регион, то существуют участки, принадлежащие отдельным потокам и недоступные всем остальным потокам. Они образуют класс частной (private) памяти потока. Частными являются:
•	стек потока, доступ к элементам которого осуществляется неявным образом при вызове функций или возврате из них;
•	переменные, объявленные внутри структурного блока, исполняемого в качестве параллельного региона (за исключением тех, которые объявлены как static).
К частной памяти относятся также регистры и локальный кэш процессора, на котором исполняется поток.
Для каждого параллельного региона могут быть указаны общие переменные, которые при создании региона преобразуются в частные. Каждая такая переменная «размножается» путем создания ее копии в частной памяти каждого потока в момент входа программы в параллельный регион. Для этих переменных могут быть указаны способы формирования начальных значений, а также возможность наследования главным потоком тех значений, которые они получили к моменту завершения параллельного региона. Могут быть также явно указаны списки общих переменных, которые используются как общие потоками внутри региона.
Особо отметим, что куча (heap) программы является общедоступной, но указатели на ее блоки могут быть именно как общими, так и частными.
Любая реализация технологии ОрепМР включает в себя следующие компоненты.
•	Директивы компилятора — используются для создания потоков, распределения работы между потоками и синхронизации потоков. Директивы пишутся как специального вида комментарии в тексте программы.
•	Библиотека функций — используется для получения и, возможно, модификации некоторых атрибутов потоков, а также для их синхронизации во время выполнения программы. Вызовы функций могут включаться программистом в текст программы.
•	Переменные окружения — используются администратором системы для управления поведением параллельной программы.
Синтаксис директив компилятора и вызова функций библиотеки времени выполнения подчиняется правилам, которые
31
различаются для разных языков программирования. Совокупность таких правил для одного языка программирования называется привязкой к языку. Технология ОрепМР создана и поддерживается для языков С, С + + и Fortran. В данном пособии рассматриваются только привязки к языкам C/C++.
Не все компиляторы поддерживают технологию ОрепМР. Если параллельную программу транслировать компилятором, не поддерживающим эту технологию, то результат трансляции будет построен в обычном последовательном (однопоточном) варианте. Для того чтобы компилятор нужным образом реагировал на директивы ОрепМР и строил параллельную программу, в командной строке его запуска должен быть указан специальный ключ. При компиляции для платформы Windows с использованием Microsoft Visual Studio включить поддержку ОрепМР можно, установив флажок «ОрепМР support» в свойствах проекта на закладке «Configuration properties\C/C++\ Language». Если приложение строится компилятором gcc для платформы Linux, то в командной строке нужно указать ключ «-fopenmp».
Чтобы в программе на C/C++ стали доступны функции библиотеки ОрепМР, в нее нужно включить заголовочный файл omp.h:
#include <omp.h>
Формат записи директив ОрепМР в программах на C/C++:
#pragma omp <имя директивы> [ опция_1[ опция_2,..]]
Без слова отр все директивы, имена которых определены в ОрепМР, смысла не имеют и будут проигнорированы компилятором. Каждая директива вместе со всеми ее опциями обязательно должна занимать ровно одну строку текста. Переносы на две или более строк не допускаются. Имеются такие директивы, которые явно не связаны с предыдущим или последующим оператором программы. Некоторые директивы определяют способ компиляции и выполнения одного оператора (или блока операторов, заключенных в фигурные скобки), непосредственно следующего за директивой в тексте программы. Он обозначается в пособии так:
<структурный блок>
32
В промежутке между директивой и ее структурным блоком могут быть записаны комментарии.
3.2.	Директивы управления регионами и потоками
3.2.1.	Директива создания параллельного региона:
#	pragma omp parallel [опция_1[ опция_2...]]
<структурный блок>
С помощью опций этой директивы можно указать:
•	требуемое количество потоков и (num_threads(n));
•	условие, при истинности которого параллельный регион действительно создается (if(logical_expression )), например if(arraySize > 100000);
•	список общих переменных для всех потоков данного региона (shared(variablesjist)), например shared(х, alpa, zeroCount);
•	список общих переменных, которые должны стать частными в каждом потоке региона и получить неопределенные начальные значения (private (variables _list)), например private (к);
•	список переменных, которые должны стать частными в каждом потоке региона и унаследовать значения из главного потока (firstprivate(variables_list)), например firstprivate(sum, ind);
•	способ назначения класса памяти по умолчанию (default (shared | попе)) всем переменным, используемым в регионе; слово попе означает, что класс памяти всех не локальных переменных (т. е. не объявленных внутри структурного блока) должен быть задан явно опциями shared, private или firstprivate; опция default(shared) может быть указана, но смысла это не имеет;
•	список глобальных (объявленных вне всех функций) переменных программы (copyin(variables_list)), перечисленных в директиве threadprivate (см. подп. 3.3.5), значения которых сохраняются в момент выхода из одного параллельного региона и должны восстанавливаться в момент входа в данный регион;
•	оператор сведения данных и список общих переменных reduction(operator : variablesjist); для каждой указанной в списке переменной создаются локальные копии в каждом потоке; локальные копии инициализируются соответственно типу
33
оператора (для аддитивных операций — ноль или его аналоги, для мультипликативных операций — единица или ее аналоги); над всеми локальными копиями каждой переменной после завершения параллельного региона будет выполнен заданный оператор сведения, результаты сведения будут занесены в одноименные общие переменные; в качестве оператора можно указывать:	|,	, &&, | | (начиная с версии 3.1
ОрепМР доступны операторы max и min); эта опция может быть очень полезна для предотвращения гонок данных; например reduction (+ : count, sum) или reduction (& : v).
Опции в любой директиве не обязательны (на это указывают квадратные скобки), если они записываются, то в произвольной последовательности и любом количестве (может быть записано несколько опций с одним и тем же именем), но смысл их не должен быть противоречивым. Это значит, что нельзя одну и ту же переменную указать, например, как общую и как частную.
Приведем пример фрагмента ОрепМР-программы, в котором используется сведение данных для параллельного вычисления суммы элементов массива. Заметим, что этот пример иллюстрирует также способ распределения итераций цикла между потоками, который может применяться и тогда, когда количество итераций заранее неизвестно.
double sum, *array;
... // формирование значений элементов массива array #pragma omp parallel reduction( + : sum)
{
int size = omp get num threads( );
int index = omp get thread num( );
while( index < count )
{
sum += array[ index ];
index += size;
}
}
3.2.2.	Директива, требующая исполнения охваченного ею участка кода в точности одним потоком:
#pragma omp single [опция_1[ опция_2...]]
<структурный блок>
Какой именно поток региона будет выполнять охваченный директивой структурный блок — не определено. Все потоки, кроме одного, обходят (не выполняют) этот фрагмент программы.
34
Опции директивы позволяют указать:
•	отмену ожидания (nowait) потоками, достигшими точки, следующей за блоком операторов, охваченным директивой single; при отсутствии этой опции все потоки региона, достигшие этой точки, останавливаются и ждут, пока не будет выполнен структурный блок;
•	списки переменных private nfirstprivate, имеющие такой же синтаксис и семантику, как в директиве parallel;
•	список переменных (copyprivate(variables_list)), значения которых после завершения выполнения структурного блока, заданного в директиве single, будут размножены во все одноименные переменные всех потоков региона; эта опция не может быть указана совместно с опцией nowait.
3.2.3.	Директива, требующая исполнения охваченного ею участка кода главным потоком программы:
#pragma omp master
<структурный блок>
Этот структурный блок будет выполнен только главным потоком программы. Остальные потоки просто пропускают данный участок и продолжают выполнение программы с оператора, расположенного следом за ним, не ожидая завершения выполнения блока мастер-потоком.
3.3.	Директивы распределения вычислений
3.3.1.	Директива распределения итераций цикла между потоками параллельного региона:
#	pragma omp for [опция_1[ опция_2...]]
<оператор цикла или гнездо операторов цикла>
В этой форме директива может быть записана внутри структурного блока директивы parallel. Директивы parallel и for можно объединять в одну директиву, с указанием совместных опций:
#	pragma omp parallel for [опция_1[ опция_2...]]
Заголовок цикла, итерации которого распределяются между потоками параллельного региона, должен выглядеть так:
35
for(var = const_el; van rel_op const_e2;
van modif_op const_e3)
и удовлетворять следующим требованиям:
•	единственной индексной переменной var должно присваиваться константное начальное значение constjel, которое может быть вычислено компилятором или исполняющей системой к моменту входа в цикл (это справедливо и для const_e2 и const_e3);
•	в качестве rel_op могут использоваться знаки >, <, >=, <=;
•	в качестве modif_op могут использоваться =, +=, -=, ++, — (две последние могут быть как в префиксной, так и в постфиксной форме, но при их использовании должно отсутствовать выражение const_e3); если используется знак простого присваивания =, то справа от него могут быть только var + const_e3, var - const_e3 или const_e3 + var.
Среди допустимых опций этой директивы есть уже известные private, firstprivate и reduction. Кроме них в директиве for могут быть указаны:
— способ распределения итераций цикла между потоками параллельной секции (schedule(type[, chunk])); параметр chunk этой опции определяет размер блока потоков, параметр type указывает тип распределения и может иметь значения:
•	static — статический (блочный, циклический или блочноциклический в зависимости от значения параметра chunk);
•	dynamic — динамический блоками размера chunk;
•	guided — динамический с экспоненциальным уменьшением размера блока от начального, формируемого автоматически, до заданного параметром chunk;
•	auto — способ распределения выбирается компилятором или исполняющей системой ОрепМР;
•	runtime — способ распределения задается значением переменной окружения OMP_SHEDULE, устанавливаемой администратором системы.
—	список переменных (lastprivate(variables_list)) главного потока, которым будут присвоены значения, вычисленные при выполнении последней итерации цикла независимо от того, какой поток эту итерацию выполнял;
—	возможность (опция ordered) появления в теле цикла директивы ordered, которая требует выполнения охваченного ею блока операторов точно в том же порядке, в котором он выполняется в последовательной версии данного цикла;
36
—	отмену (nowait) ожидания всеми потоками полного завершения выполнения всех итераций цикла; без этой опции никакой поток не может выйти за пределы цикла, пока не будут выполнены все итерации; эта опция не может быть указана совместно с опцией lastprivate;
—	глубину п вложенных друг в друга циклов гнезда (collapse(n)), пространство итераций которых подлежит распределению между параллельными потоками (доступно только в реализациях версий 2.5 и более поздних).
Перепишем фрагмент программы, вычисляющий сумму элементов массива из подп. 3.2.1 с использованием директивы for:
double sum, *array;
int i;
... // формирование значений элементов массива array
#	pragma omp parallel for reduction( + : sum)
for(i = 0; i < count; i++) {
sum += array[ i ];
}
3.3.2.	Директива, указывающая на необходимость исполнения охватываемого ею участка кода в том порядке, в котором он выполняется в последовательной программе:
#pragma omp ordered
<структурный блок>
3.3.3.	Директива задания области нециклического параллелизма (участки структурного блока кода, которые могут выполняться параллельно, выделяются с помощью описываемой далее директивы section):
#pragma omp sections [опция_1[ опция_2___]]
<структурный блок>
В качестве опций этой директивы можно указывать private, firstprivate, lastprivate, reduction и nowait, имеющие в точности такой же синтаксис и тот же смысл, что и у директивы fo г.
3.3.4.	Директива определения участка нециклического кода для одного потока:
#pragma omp section
<структурный блок>
37
С помощью этой директивы внутри участка кода, охваченного директивой sections, определяются отдельные фрагменты для исполнения параллельными потоками. Перед самым первым фрагментом участка нециклического кода после директивы sections директиву section можно не писать.
Связанная пара директив sections и section может быть использована, когда в программе можно выделить различные информационно независимые (или слабо зависимые) времяза-тратные фрагменты.
Например, пусть есть «поставщик» блоков данных и их «потребитель». Алгоритмы работы поставщика и потребителя существенно различаются, при этом и тот и другой в определенной степени зависят от внешнего оборудования, т. е. могут инициироваться или приостанавливаться некоторыми внешними событиями. В этих условиях можно применить распараллеливание с использованием директив sections и section и структуры типа vector для передачи указателей на блоки данных между поставщиком и потребителем:
vector<dataBlock*> queue;
int pushed = 0;
int pulled = 0;
#pragma omp parallel sections
{
#pragma omp section
while(true) { // это поставщик
dataBlock* dBlockPtr;
___ //формирование блока данных и указателя dBlockPtr queue.push(dBlockPtr);
pushed += 1;
}
#pragma omp section
while(true) { // это потребитель
if( pushed > pulled) {
dataBlock* dBIPtr = queue[pulled];
pulled += 1;
//обработка блока данных по указателю dBIPtr
}
}
}
Нужно отметить, что в этом фрагменте отсутствуют действия по удалению обработанных указателей из вектора queue.
38
В связи с тем, что векторы не являются потокобезопасными структурами данных, реализация удаления элементов из них в параллельной программе не является простой задачей. Для ее решения необходимо обеспечивать синхронизацию потоков, реализующих алгоритмы поставщика и потребителя. Средства синхронизации рассматриваются в пп. 3.4 и 3.6.
3.3.5.	Директива объявления списка локальных переменных потоков, сохраняемых между параллельными регионами:
#pragma omp threadprivate( variables_list )
Эта директива позволяет превратить перечисленные в списке variables_list статические переменные (которые по умолчанию являются общими) в частные переменные, значения которых сохраняются при закрытии каждого параллельного региона и восстанавливаются для каждого потока в момент открытия следующего.
3.3.6.	Директива создания отдельной независимой задачи (появилась в стандарте версии 2.5 ОрепМР):
#	pragma omp task [опция_1[ опция_2...]]
<структурный блок>
Этой директивой поток оформляет указанный структурный блок в качестве отдельной задачи, которая ставится системой поддержки ОрепМР в очередь задач и будет выполнена позже каким-либо, возможно тем же самым, потоком. Поток, создавший задачу, «перешагивает» через структурный блок, но позже может остановиться до завершения всех созданных им задач (см. подп. 3.3.7).
Возможные опции:
•	if, default, private, firstprivate и shared — имеют такой же синтаксис и смысл, как и в предыдущих директивах;
•	untied — означает, что задача может быть выполнена любым потоком данного региона; если опция не указана, то задача может быть выполнена только породившим ее потоком (это поведение может быть иным и зависит от реализации).
3.3.7.	Директива ожидания потоком завершения всех независимых задач, созданных именно данным потоком:
#pragma omp taskwait
39
3.3.8.	Директива возможной приостановки выполнения задачи для того, чтобы дать возможность выполнения других задач:
#pragma omp taskyield
Директивы task, taskwait, и taskyield являются наиболее удобным средством программирования рекурсивного (неитеративного) параллелизма. Например, пусть нужно реализовать параллельный обход двоичного дерева. Это значит, что начиная с корня дерева нужно обойти его левое поддерево, затем правое поддерево и, наконец, обработать содержимое корня. Для каждого поддерева должны быть выполнены те же самые действия, поэтому процесс обхода носит выраженный рекурсивный характер, который невозможно распараллелить с помощью директивы for (в том числе и потому, что количество итераций этого циклического процесса заранее неизвестно).
Задачи обхода деревьев встречаются довольно часто, поэтому рассмотрим возможный способ их решения.
Пусть узел дерева определен так:
struct node {
struct node* leftSubTree;
struct node* rightSubTree;
struct content { ... } content;
};
Параллельный обход дерева с корнем treeRoot может быть реализован таким фрагментом:
#pragma omp parallel
{
tfpragma omp single
{
#	pragma omp task firstprivate( treeRoot ) untied subTreeTraversal( treeRoot );
tfpragma omp taskwait
}
}
Создаваемая любым (одним) потоком параллельного региона задача состоит в вызове функции subTreeTraversal. Эта функция реализует создание новых задач для рекурсивного обхода поддеревьев переданного ей узла (если они есть), дожидается 40
их завершения и вызывает функцию processNode, выполняющую обработку содержимого корня дерева/поддерева:
void subTreeTraversal(struct node* root) {
if(root -> leftSubTree)
#	pragma omp task firstprivate(root -> leftSubTree) untied subTreeTraversal(root -> leftSubTree);
if(root -> leftSubTree)
#	pragma omp task firstprivate(root -> rightSubTree) untied subTreeTraversal(root -> rightSubTree);
#	pragma omp taskwait processNode(root);
}
Все создаваемые задачи могут выполняться любым потоком региона. Поток, инициировавший создание самой первой задачи, дожидается завершения ее выполнения, а значит — завершения выполнения всех подзадач, созданных ею. Каждая из них, в свою очередь, ожидает завершения созданных ею подзадач и т. д.
3.4.	Директивы синхронизации
Синхронизация действий потоков необходима во всех случаях, когда они тем или иным способом модифицируют общие данные, т. е. когда может возникнуть состояние гонки данных.
3.4.1.	Директива явной барьерной синхронизации:
#pragma omp barrier
Потоки, выполняющие текущий параллельный регион, дойдя до этой директивы, переводятся системой ОрепМР в состояние ожидания. Выход из этого состояния осуществляется в тот момент, когда директиву barrier выполнит последний (по времени) поток региона.
Следует отметить, что неявно барьерная синхронизация выполняется при завершении параллельных регионов, циклов и секций в случае, если не использована опция nowait.
3.4.2.	Директива, объявляющая критическую секцию — участок кода программы, который единовременно может выполняться не более чем одним потоком:
#pragma omp critical [(< section_name >)]
<структурный блок>
41
Если критическая секция уже выполняется каким-либо потоком, то все остальные, выполняющие эту директиву для секции с данным именем, будут заблокированы до тех пор, пока находящийся в секции поток не закончит выполнение структурного блока кода. Как только это произойдет, один из заблокированных потоков войдет в критическую секцию. Если в состоянии ожидания входа в критическую секцию находилось несколько потоков, то случайным образом выбирается один из них, остальные заблокированные продолжают ожидание.
Все неименованные критические секции ассоциируются с одним и тем же пустым именем. Все критические секции, имеющие одно и то же имя, рассматриваются как одна секция, даже если находятся в разных параллельных областях. Побочные входы в структурный блок и выходы из него запрещены.
Критические секции с одним и тем же именем обычно используются для блокировки возможного одновременного доступа нескольких потоков к одному общему ресурсу (переменной, структуре, экземпляру класса и т. д.) с целью его модификации. Этот механизм позволяет предотвращать гонки данных, но является довольно времяемким. Во многих случаях более предпочтительным является использование директивы atomic.
3.4.3.	Директива блокировки доступа к общей переменной из левой части оператора присваивания на время выполнения всех действий с этой переменной в данном операторе:
#pragma omp atomic
<оператор присваивания>
Атомарной (блокирующей выполнение остальных потоков) является только работа с переменной из левой части оператора присваивания, при этом вычисления в его правой части, использующие другие переменные, могут выполняться одновременно в нескольких потоках.
3.4.4.	Директива актуализации значений переменных:
#pragma omp flush [(variables_list)]
Значения всех переменных из списка variables_list, если он задан (или всех переменных потока), временно хранящиеся в регистрах и кэш-памяти исполняющего поток ядра, переносятся в основную память. Следовательно, все изменения пере
42
менных, сделанные потоком до исполнения этой директивы, станут видимы остальным потокам. Эти действия производятся только с данными выполнившего директиву потока, переменные, изменявшиеся другими потоками, не затрагиваются. До полного завершения этих действий никакие другие операции с участвующими в директиве flush переменными не могут выполняться. Поэтому выполнение директивы flush без списка переменных может повлечь значительные накладные расходы. Если в данный момент нужна гарантия согласованного представления не всех, а лишь некоторых переменных, то именно их следует явно перечислить в директиве списком.
3.5.	Переменные окружения
На то, как будет исполняться параллельная программа, могут влиять не только директивы ОрепМР, но и переменные окружения, значения которых могут устанавливаться или изменяться администратором системы. Некоторые переменные могут модифицироваться самой параллельной программой путем вызова соответствующих функций библиотеки (см. п. 3.6). В список переменных окружения входят следующие переменные.
OMP_NUM_THREADS — определяет количество потоков в новой параллельной области.
OMP_NESTED — определяет возможность создания вложенных параллельных регионов.
OMP_MAX_ACTIVE_LEVELS — задает максимальную глубину вложенности параллельных регионов.
OMP_DYNAMIC — определяет возможность динамической установки количества потоков в регионе (имеет больший приоритет, чем OMP_NUM_THREADS).
OMP_THREAD_LIMIT — задает максимальное количество потоков для всей программы.
OMP_SCHEDULE — содержит способ распределения итераций цикла (используется, если указано значение runtime опции schedule директивы for).
OMP_STACKSIZE — устанавливает размер стека каждого потока ОрепМР-программы.
OMP_WAIT_POLICY — разрешает или запрещает выделять кванты процессорного времени ждущим потокам.
43
З.б.	Библиотека функций ОрепМР
Функции библиотеки позволяют получать и изменять некоторые параметры регионов и потоков (обычно являющиеся значениями переменных окружения).
double omp_get_wtime(void); — возвращает текущее время.
double ompjget_wtick(void); — возвращает размер тика таймера.
int omp_get_thread_num(void); — возвращает номер потока.
int ompjgetjiumjdireads (void); — возвращает количество потоков в текущем регионе (1, если вызывается вне региона).
int omp_get_max_threads(void); — возвращает максимально возможное количество потоков для следующей параллельной области.
int omp_get_thread_limit(void); — возвращает значение ОМР_ THREAD_LIMIT.
void omp_set_num_threads(int num); — устанавливает новое значение переменной OMP_NUM_THREADS.
int omp_get_numj>rocs(void); — возвращает количество про-цессоров/ядер.
int omp_get_dynamic(void); — возвращает значение ОМР_ DYNAMIC.
void omp_set_dynamic(int num); — устанавливает новое значение переменной OMP_DYNAMIC.
int omp_get_nested(void); — возвращает значение переменной OMP_NESTED.
void omp_set_nested(int nested); — устанавливает новое значение переменной OMP_NESTED.
int omp_in_parallel(void); — возвращает 0, если функция вызвана из последовательной области, и 1, если из параллельной.
int omp_get_max_active_levels(void); — возвращает значение переменной OMP_MAX_ACTIVE_LEVELS.
void omp_set_max_active_levels(int max); — устанавливает новое значение OMP_MAX_ACTIVE_LEVELS.
int omp_get_level(void); — возвращает глубину вложенности параллельных областей.
int omp_get_ancestor_thread_num(int level); — возвращает номер потока, породившего текущую параллельную область.
int omp_get_team_size(int level); — возвращает для заданного параметром level уровня вложенности параллельных областей количество потоков, порожденных одним родительским потоком.
44
int omp_get_active_level(void); — возвращает количество вложенных параллельных областей, обрабатываемых более чем одним потоком.
3.6.1.	Следующая группа функций используется для синхронизации доступа параллельно выполняющихся потоков к общим ресурсам с помощью так называемых замков — общих переменных типа omp_lock_t или omp_nest_lock_t.
void omp_init_lock(pmp_lock_t ★lock); — создать (инициализировать) простой замок.
void omp_init_nest_lock(pmp_nest_lock_t *lock); — создать замок с множественными захватами.
void omp_destroy_lock(pmp_lock_t *lock); — уничтожить простой замок (перевести в неинициализированное состояние).
void omp_destroy_nest_lock(pmp_nest_lock_t *lock); — уничтожить множественный замок (перевести в неинициализированное состояние).
void omp_set_lock(pmp_lock_t ★lock); — захватить простой замок (если замок уже захвачен другим потоком, то данный поток переводится в ждущее состояние).
void omp_set_nest_lock(pmp_nest_lock_t ★lock); — захватить множественный замок (если замок уже захвачен другим потоком, то данный поток переводится в ждущее состояние; если замок захвачен этим же потоком, то увеличивается счетчик захватов).
void omp_unset_lock(pmp_lock_t ★lock); — освободить простой замок (если есть потоки, ждущие его освобождения, то один из них, выбираемый случайным образом, захватывает этот замок и переходит в состояние выполнения).
void omp_unset_nest_lock(omp_lock_t ★lock); — уменьшить на 1 количество захватов множественного замка (если количество захватов стало равно нулю и есть потоки, ждущие его освобождения, то один из них, выбираемый случайным образом, захватывает этот замок и переходит в состояние выполнения).
int omp_test_lock(pmp_lock_t ★lock); — попытаться захватить простой замок (если попытка не удалась, т. е. замок уже захвачен другим потоком, возвращается 0, иначе возвращается 1).
int omp_test_nest_lock(pmp_lock_t ★lock); — попытаться захватить множественный замок (если попытка не удалась, т. е. замок уже захвачен другим потоком, возвращается 0, иначе возвращается новое значение счетчика захватов).
45
Следует помнить, что замки являются наиболее медленным средством синхронизации процессов, быстрее работают критические секции, еще быстрее — директива atomic, и, если ее можно применить — опция reduction в директивах управления параллелизмом.
Контрольные вопросы и задания
1.	Из каких составных частей состоит технология ОрепМР?
2.	С помощью какой директивы (директив) создаются новые параллельные регионы программы?
3.	Что такое критическая секция программы?
4.	Каким образом можно установить нужное количество потоков для создания очередного параллельного региона?
5.	Как обеспечить выполнение фрагмента параллельного региона только главным потоком?
6.	Какая опция директивы ОрепМР for используется для указания способа распределения итераций цикла между потоками параллельного региона?
7.	Что такое deadlock? Каким правилам нужно следовать, чтобы избежать возможности попадания параллельной программы в deadlock?
8.	Что такое сведение данных? Какие опции и в каких директивах используются для выполнения сведения данных?
9.	Что делает директива ОрепМР threadprivate?
10.	Как обеспечить выполнение фрагмента параллельного региона потоком с максимальным номером в данном параллельном регионе?
11.	Для чего используется опция firstprivate? Чем она отличается от опций private и lastprivate?
12.	Что такое вложенный параллельный регион программы? В каких случаях его нельзя создать?
13.	Что такое неявная барьерная синхронизация? С помощью каких средств ее можно отменить?
14.	Для чего используются директивы ОрепМР sections и section? Что делает каждая из этих директив?
15.	Перечислите все средства синхронизации потоков в ОрепМР.
16.	Перечислите возможные способы распределения итераций цикла между потоками.
17.	Что делает директива ОрепМР atomic? Можно ли одну и ту же переменную записать и в левой и в правой частях оператора присваивания после этой директивы?
18.	В чем состоит различие между общими и локальными переменными потока?
19.	С помощью каких средств можно ограничить глубину вложенности параллельных регионов программы?
20.	От чего зависит равномерность загрузки процессоров/ядер системы с общей памятью?
46
21.	Каким образом функция, вызываемая из параллельной программы, может выяснить, в последовательном или параллельном регионе она выполняется?
22.	Как обеспечить выполнение фрагмента параллельного региона в точности одним потоком?
23.	Могут ли два потока одновременно захватить один множественный замок?
24.	Для чего может использоваться опция reduction в некоторых директивах ОрепМР?
25.	Каково назначение директивы ОрепМР single?
26.	Каким ограничениям должен удовлетворять цикл, для того чтобы к нему можно было применить директиву ОрепМР for?
27.	Чем отличается простой замок от множественного?
28.	Перечислите различия механизмов синхронизации потоков в ОрепМР.
29.	Перечислите способы распределения итераций цикла между потоками параллельного региона.
Тема 4
CUDA: ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ НА ГРАФИЧЕСКОМ ПРОЦЕССОРЕ
CUDA (Compute Unified Device Architecture) представляет собой совокупность аппаратного и программного обеспечения компании NVIDIA, обеспечивающая возможность подготовки и исполнения высокопараллельных программ на графических процессорах [14—17].
4.1.	Архитектура CUDA
Графический процессор (Graphics Processor Unit, GPU) — это специализированное вычислительное устройство, которое является сопроцессором к основному процессору компьютера (central processor unit, CPU). Графический процессор реализует обработку данных согласно модели ОКМД, имеет собственную иерархически организованную память и способен параллельно выполнять очень большое количество (тысячи и десятки тысяч) отдельных потоков управления.
Технология ОрепМР, рассмотренная в теме 3, и технология CUDA вполне совместимы и могут использоваться в одной параллельной программе для достижения максимально возможного ускорения вычислений. При этом между потоками, выполняемыми на ядрах CPU, и потоками, выполняемыми графическим процессором, есть принципиальные различия:
•	потоки графического процессора обладают крайне низкой «стоимостью»: создание и управление ими требует значительно меньших ресурсов центрального процессора по сравнению с потоками ОрепМР;
•	для эффективной утилизации возможностей GPU нужно запускать многие тысячи и даже миллионы потоков, тогда как для ядер CPU обычно бывает трудно организовывать больше нескольких десятков потоков.
48
Приложения, использующие возможности CUDA для параллельной обработки данных, взаимодействуют с GPU через программные интерфейсы, называемые CUDA-runtime и CUDA-driver (обычно CUDA-runtime и CUDA-driver одновременно не используются, хотя в принципе это возможно).
CUDA-программы пишутся на «расширенном» языке C/C++, при этом их параллельная часть (это так называемые функции-ядра) выполняется на графическом процессоре, а обычная последовательная часть — на центральном процессоре. Компоненты архитектуры CUDA, принимающие участие в подготовке и исполнении приложения, автоматически осуществляют разделение частей и управление их запуском.
4.2.	Аппаратные компоненты CUDA
Графический процессор NVIDIA состоит из нескольких так называемых кластеров текстурных процессоров (texture processor cluster), как показано на рис. 4.1.
49
Каждый такой кластер состоит из логики работы с текстурами ЛТ (она в параллельных вычислениях не участвует и далее не рассматривается) и двух или более потоковых мультипроцессоров ПМ (Streaming Multiprocessor, SM). Каждый из них содержит:
•	начало конвейера (front end), выполняющее чтение и декодирование инструкций, а также передачу их на выполнение; элементами начала конвейера являются кэш команд, кэш данных и планировщик потоков (блок выборки инструкций);
•	конец конвейера (back end), состоящий из нескольких (от восьми и более) потоковых процессоров ПП (Streaming Processor, SP) и специализированных функциональных устройств СФУ (Special Function Unit, SFU).
Потоковые процессоры «умеют» выполнять обычные арифметические операции над целыми числами и числами с плавающей точкой одинарной точности (32 бита). Специализированные функциональные устройства выполняют быстрые приближенные трансцендентные операции, такие как вычисление встроенных функций___cosf(), _expf() и многих других.
Минимальной единицей работы для потокового мультипроцессора является так называемый ворп (warp) — пучок из 32 потоков. Потоки пучка идентифицируются индексами, увеличивающимися на единицу, т. е. являются близкими «соседями». Пучки группируются в блоки таким образом, чтобы на выполнение блока хватило ресурсов одного потокового мультипроцессора. Пучки могут быть готовыми к выполнению или находиться в состоянии ожидания завершения некоторой запущенной ранее операции (например — конкурентного обращения к памяти). Если хотя бы один поток ворпа чего-то ждет, то весь ворп находится в состоянии ожидания.
В потоковом мультипроцессоре инструкции выполняются по принципу SIMD, то есть одна инструкция применяется ко всем потокам в ворпе, но NVIDIA называет такой способ выполнения SIMT (single instruction multiple threads, одна инструкция, много потоков). Важно отметить, что конец конвейера работает на частоте, в два раза превосходящей частоту работы начала конвейера. На практике это означает, что данная часть выглядит в два раза «шире», чем она есть на самом деле (то есть как 16-процессорная SIMD-совокупность вместо восьмипроцессорной) .
50
На каждом такте работы (такте своей частоты) планировщик потоков выбирает ворп, готовый к выполнению, и запускает выполнение очередной его инструкции. Некоторые инструкции (например, сложение, умножение и т. п.) выполняются потоковыми процессорами, некоторые — специализированными функциональными устройствами. Поскольку те и другие могут работать одновременно, наилучшим с точки зрения максимума производительности и загрузки оборудования является чередование обычных операций и вызовов встроенных функций.
Однако такой полностью параллельный режим работы потокового мультипроцессора может нарушаться и при выполнении операций ветвления, которые могут присутствовать в коде исполняемой программы. После выполнения любой операции ветвления планировщик приостанавливает все потоки одной ветки и организует выполнение потоков другой ветки до их перехода в состояние ожидания или до завершения. После этого возобновляются приостановленные потоки. В силу этого то одна, то другая часть потоковых процессоров может простаивать.
Ворпы являются составной часть более крупного конгломерата потоков — блока. Объединение ворпов в блоки оказалось целесообразным в силу того, что тактовая частота основной памяти графического процессора значительно ниже тактовой частоты потоковых процессоров. Поэтому при обращении хотя бы одного потока какого-либо ворпа к этой памяти можно занять мультипроцессор выполнением другого готового ворпа. Тем самым реализуется «скрытие» временных затрат на доступ к медленной основной памяти. В связи с тем, что каждому потоку требуется некоторое, пусть небольшое, количество ресурсов мультипроцессора, которые должны сохраняться в течение периодов ожидания, максимальное количество потоков в блоке ограничено сверху. Для большинства графических процессоров NVIDIA количество потоков в блоке не должно превышать 1024 (для младших моделей — 512). Блок потоков всегда исполняется одним мультипроцессором от момента начала до момента завершения всех потоков, т. е. не может перемещаться во время выполнении ядра на другой мультипроцессор. Размер блока может устанавливаться при запуске параллельной части в диапазоне от 1 до 1024 потоков и должен выбираться программистом в зависимости от множества параметров, в том числе — от необходимого каждому потоку количества внутрен
51
них ресурсов мультипроцессора. К таким ресурсам, прежде всего, относится память.
Вся совокупность вычислительных процессов, запускаемых на графическом процессоре для исполнения параллельной части программы в форме функции-ядра, называется сеткой блоков потоков.
Очень важным для эффективного программирования графического процессора является понимание структуры его иерархически организованной памяти и параметров ее составных частей. Возможности потоков каждого блока, исполняемые потоковыми процессорами, по доступу к разным уровням памяти показаны на рис. 4.2. Стрелки внизу обозначают возможность доступа хост-компьютера к глобальной памяти и памяти констант как по чтению, так и по записи. Все прочие уровни памяти графического процессора для хост-компьютера недоступны.
Рис. 4.2. Структура памяти графического процессора
Глобальная память доступна любому потоку как на чтение, так и на запись, имеет объем от нескольких сотен мегабайт
52
до десятков гигабайт и относительно большое время доступа. Поэтому в каждом потоковом мультипроцессоре обычно имеется кэш памяти, не показанный на рис. 4.2.
Каждый потоковый мультипроцессор (на рис. 4.2 показанный как блок потоков) имеет общую для всех потоков блока разделяемую память объемом 32—48 килобайт. Это значительно более быстрая память, чем глобальная, время обращения к ней в 300—500 раз меньше. Разделяемая память каждого блока недоступна для потоков любого другого блока. При программировании графического процессора рекомендуется учитывать различие скоростей разных уровней памяти и перемещать переменные, к которым часто обращаются потоки, из главной в разделяемую память для достижения максимальной производительности.
Потоковый мультипроцессор имеет также так называемый регистровый файл — быструю память объемом 64 килобайта, которая поровну распределяется между потоками блока. В этой памяти хранятся частные данные потока, поэтому доступа к ней не имеет никакой другой поток. Требуемый потоку объем частной памяти определяется тем, сколько локальных переменных объявлено программистом в тексте функции-ядра. Если этот объем больше, чем размер регистрового файла, деленный на количество потоков в блоке, то CUDA обеспечивает захват недостающего количества байт на каждый поток в глобальной памяти и размещает в ней некоторые частные переменные потока. Такие участки глобальной памяти, доступные теперь только потокам-владельцам, называются локальной памятью (ЛП на рис. 4.2). Доступ к этим частным переменным будет существенно медленнее, чем к переменным, размещенным в регистрах. Об этом следует помнить и при запуске ядра заботиться о правильной установке количества потоков в блоке для достижения максимальной производительности.
Архитектура CUDA предусматривает широкий набор возможностей и поддерживается в видеокартах нескольких линеек (GEforce, Quadro, NVS, Tesla, Tegra, Jetson). Разные модели графических процессоров, функционирующих в таких видеокартах, имеют различные наборы параметров (отличающиеся количеством регистров, объемом разделяемой памяти мультипроцессора) и поддерживают разные наборы возможностей.
53
4.3.	Программные компоненты CUDA
Программная составляющая CUDA обеспечивает:
•	компиляцию исходных текстов программ с формированием последовательной и параллельной частей приложения;
•	управление графическим процессором по запросам из последовательной части приложения для выполнения на нем параллельной части;
•	подключение библиотек функций параллельных расчетов;
•	профилирование и отладку параллельной части приложения.
Для программирования графических процессоров специалистами фирмы NVIDIA был синтаксически расширен базовый язык C/C++. Преобразование исходных текстов, написанных на расширенном варианте языка и содержащих как последовательную, так и параллельную части, осуществляет компилятор nvcc, использующий при этом компилятор и компоновщик базовой платформы для формирования кода приложения. Схема процесса компиляции показана на рис. 4.3. По существу, nvcc выполняет препроцессинг, выделяя и обрабатывая только параллельную часть приложения и передавая последовательную часть на обработку компилятору базовой платформы (например, gcc в Linux).
После запуска сформированного компоновщиком приложения оно взаимодействует с аппаратурой графического процессора, вызывая либо функции CUDA-драйвера, либо функции CUDA-рантайма, либо функции специализированных библиотек, как показано на рис. 4.4. Все функции CUDA-драйвера имеют префикс «си» (например, cuMemAlloc, cuLaunchKernel и т. д.), все функции CUDA-рантайма — префикс «cuda» (например, cudaGetDeviceProperties, CudaMemCopy и т. д.). Функции CUDA-библиотек специфичны и общего классификатора не имеют.
Профилирование и отладка CUDA-приложения могут выполняться обычными средствами, но только до тех пор, пока исследуются процессы последовательной части программы. Для профилирования и отладки параллельной части компания NVIDIA разработала инструментарий NVIDIA Parallel Nsight, интегрирующийся с Microsoft Visual Studio. Он работает либо на хост-компьютере с двумя графическими процессорами (ло
54
кальный вариант), либо на двух компьютерах, каждый из которых оснащен одним графическим процессором (удаленный вариант).
Рис 43. Схема процесса компиляции CUDA-программ
В момент запуска на выполнение параллельная часть приложения (ядра) организуется как сетка (grid) блоков (block) потоков (thread). Потоки одного блока могут взаимодействовать через разделяемую память и синхронизироваться между собой. Потоки, принадлежащие разным блокам, взаимодействовать друг с другом могут только через глобальную память устройства. Однако вследствие того, что не существует возможности синхронизации работы разных блоков, организация таких взаимодействий связана с большими трудностями и их следует избегать.
Сетка (и независимо от нее блоки) могут быть организованы в одно-, двух- или трехмерную структуру блоков и пото
55
ков соответственно. Количество измерений и размеры (сетки в блоках, а блоков — в потоках) по каждому из них для каждого запуска ядра на графическом процессоре задает программист. Каждый поток после запуска имеет доступ к переменным (структурам), хранящим его собственные координаты внутри блока и координаты содержащего поток блока внутри сетки по каждому из измерений. Это позволяет распределять вычислительные работы между потоками на основании (или с использованием) координатных индексов потока. Возможность формировать n-мерные структуры сетки блоков потоков позволяет гибко адаптировать разрабатываемое приложение к природе решаемой задачи. Например, при решении задач наподобие сортировки массивов следует запускать одномерную сетку блоков потоков, в матричной алгебре естественной будет двумерная организация сетки и блоков и т. д.
Рис. 4.4. Схема взаимодействия приложения с компонентами CUDA
4.4.	CUDA-расширение языка С/С++
Расширения языка C/C++ включают:
•	спецификаторы функций, указывающие, где должна выполняться определяемая функция и откуда она может быть вызвана;
•	спецификаторы переменных, задающие тип памяти, используемый для определяемой переменной;
56
•	дополнительные типы данных.
•	директива запуска ядра, задающая его аргументы и структуру сетки блоков потоков;
•	встроенные переменные ядра, содержащие координаты потока и блока.
4.4.1.	Дополнительные спецификаторы видов функций перечислены в табл. 4.1.
Таблица 4.1
Спецификатор	Исполняется на	Вызывается из
	host	<return_type> <name>0	CPU	CPU
	global	void <kernel_name>0	GPU	CPU
	device	<return_type> <name>0	GPU	GPU
Особенности дополнительных видов функций:
•	функция, объявленная как__global_(ядро), не может
возвращать никакого значения (тип возвращаемого значения — всегда void);
•	функция вида___device_и выполняется на GPU и вы-
зывается из функций, исполняемых на GPU, поэтому в хост-программе нельзя получить указатель на такую функцию.
В функциях, исполняемых на GPU:
•	локальные переменные не могут объявляться с использованием ключевого слова static;
•	количество аргументов не может быть переменным.
4.4.2.	Дополнительные спецификаторы видов переменных перечислены в табл. 4.2.
Таблица 4.2
		Тип памяти	Область видимости
	device		=shared=<type>	разделяемая	блок
	device		<type>	глобальная	сетка
	device		=constant=<type>	кэш констант	сетка
Ключевое слово__device_необязательно, если используется
__shared_или____constant_. Локальные переменные функций, выполняемых на GPU, объявляются без идентификатора вида и хранятся в регистрах. Исключения — большие структуры или массивы, размещаемые компилятором в медленной локальной памяти, поэтому их использования следует избегать.
57
4.4.3.	Дополнительные встроенные векторные типы:
[u]char[l..4], [u]short[l..4], [u]int[l..4],
[u]long[l. .4], float[1..4]
Буква и в этих объявлениях необязательна, она является сокращением от слова unsigned и в случае ее указания означает, что все поля данного типа являются беззнаковыми целыми.
Каждый такой тип является полным аналогом структуры, например, объявление
ulong4 van;
эквивалентно последовательности операторов вида:
struct unsignedLong4 {
unsigned long x;
unsigned long y;
unsigned long z;
unsigned long w;
};
struct unsignedLong4 var;
Доступ к полям переменных, имеющих эти дополнительные типы, осуществляется по именам х, у, z, w, например:
uint4 param;
int у = param.у;
unsigned int abc = param.w;
Тип dim3 является полным синонимом uint3. Обычно переменные типа dim3 используются для задания параметров сетки блоков потоков при запуске ядра.
4.4.4.	Ядро (собственно программа для GPU) и его запуск.
Прототип функции, являющейся ядром, выглядит так:
global void <kernel_name>(<formal_arguments_list>);
В хост-программе ядро запускается на выполнение с помощью специальной синтаксической конструкции, входящей в расширение языка C/C++:
kernel_name«<grid , block [, mem[, stream] ]»> ( params );
58
Здесь:
•	kernel_name — имя функции-ядра, т. е. объявленной со спецификатором____global_;
•	grid — размерность сетки в блоках, задаваемая с помощью значения типа dim3 (для одномерной сетки можно указать скалярное значение, например 256, оно автоматически будет преобразовано к типу dim3{256,1,1});
•	block — размерность блока в потоках, тоже задаваемая с помощью значения типа dim3 (также допускается возможность указания скалярного значения);
•	тет — количество дополнительной разделяемой памяти в байтах для каждого блока (выделяется в дополнение к статически объявленным разделяемым переменным в ядре и отображается на переменные, объявляемые как extern__shared___...);
этот параметр может отсутствовать, но тогда, если указан следующий параметр (stream), то после параметра block должны следовать две запятые;
•	stream — идентификатор (handle) очереди работ CUDA-stream, в котором будет работать это ядро; если этот параметр не указан, то ядро ставится в очередь работ по умолчанию, создаваемую автоматически; очереди работ не по умолчанию могут создаваться путем вызова соответствующих функций;
•	params — значения фактических аргументов функции-ядра.
4.4.5.	Встроенные переменные, создаваемые и инициализируемые автоматически при запуске ядра:
dim3 gridDinv, — размеры сетки в блоках;
dim3 blockDim; — размеры одного блока в потоках;
dim3 blockldx; — индекс блока внутри сетки;
dim3 threadldx; — индекс потока внутри блока.
Эти переменные доступны только для кода параллельных ветвей ядра, исполняемого на GPU.
Приведем простейший пример оформления функции-ядра и его запуска для вычисления поэлементной суммы двух массивов. В последовательной программе для этого необходим цикл, перебирающий все элементы обрабатываемых массивов (первая строка фрагмента — определение массивов):
float arrayOne[size], arrayTwofsize], resultfsize];
... // формирование значений массивов arrayOne и arrayTwo
for(i = 0; i < size; i += 1)
result[i] = arrayOne[i] + arrayTwo[i];
59
При преобразовании этого алгоритма для вычислений на графическом процессоре (пока будут опущены некоторые детали работы с памятью) тело цикла без его заголовка превращается в такое тело функции-ядра (параллельная часть программы):
__global__ void vAdd(int n, float *al, float *a2, float *a3) { unsigned int tid = blockDim.x * blockldx.x + threadldx.x; if(tid < n)
a3[tid] = al[tid] + a2[tid];
}
В хост-программе после формирования элементов суммируемых массивов и переноса их в память графического процессора ядро запускается так:
unsigned int grid = size / 512 + (size % 512 == 0? 0 : 1);
vectorAdd <<< grid, 512 >>> (
size, arrayOne, arrayTwo, result);
В этой последовательности операторов указано, что должна быть запущена одномерная сетка из size /512 (или size / 512 + + 1) блоков, каждый из которых тоже одномерен и состоит из 512 потоков. Одномерность определяется тем, что вместо значений типа dim3 в операторе запуска ядра указаны скаляры.
Если размер массива не делится нацело на размер блока потоков, то потоки с индексами, превышающими размер массива, не должны обращаться к памяти. Поэтому в коде ядра записано, что каждый поток прежде всего вычисляет собственную координату (индекс потока tid) по единственному измерению сетки и проверяет, находится ли этот индекс в пределах границ изменения индексов массива. Вычисление суммы двух элементов выполняется только тогда, когда это условие истинно.
4.5.	Технология разработки CUDA-программ
Технология параллельного программирования графического процессора, рекомендуемая разработчиками архитектуры CUDA и преследующая целью эффективную загрузку оборудования и достижение максимальной производительности, со
60
стоит в реализации следующих действий в последовательной и параллельных частях программы.
1.	Получить от CUDA-драйвера сведения о количестве и характеристиках графических процессоров, подключенных к хост-компьютеру; если их несколько, то нужно принять решение о том, какое именно оборудование будет использоваться.
2.	Подготовить данные для обработки и загрузить их в глобальную память GPU.
3.	Спланировать разбиение обрабатываемых данных на фрагменты, целиком помещающиеся в разделяемую память графического процессора, т. е. написать функцию-ядро.
4.	Запустить ядро. В программе ядра для каждого блока:
1)	выполнить перемещение нужных фрагментов данных из глобальной памяти в разделяемую память блока;
2)	выполнить требуемые вычисления, используя данные, хранящиеся в быстрой разделяемой памяти блока;
3)	скопировать сформированные результаты обработки из разделяемой памяти в глобальную память устройства.
5.	Перенести полученные результаты из памяти графического процессора в основную память компьютера, выполнить их окончательную обработку или вывести на внешний носитель.
6.	При необходимости повторить шаги 2—5, или 3—5, или А—5 нужное количество раз.
В качестве иллюстрации этого процесса и его результатов приведем CUDA-программу для вычисления произведения матриц. В параллельной части этой программы каждый блок потоков для вычисления собственной части матрицы результата перемещает в разделяемую память потокового мультипроцессора квадратные фрагменты полос строк матрицы-множимого (обозначены как FA1, FA2, ...) и полос столбцов матрицы-множителя (FBI, FB2, ...).
При обработке каждой пары фрагментов производится накопление частичных сумм произведений элементов умножаемых матриц, выбираемых из быстрой разделяемой памяти. Полосы строк и столбцов перемножаемых матриц, обрабатываемые некоторым блоком потоков, а также вычисляемая им подматрица результата выделены на рис. 4.5. Внутри этих полос серым цветом выделены строка и столбец, обрабатываемые одним потоком блока. Черным цветом показан элемент результирующей матрицы, формируемый этим потоком.
61
Рис 4.5. Перемножение матриц на CUDA
// Матрицы в памяти GPU хранятся в обычном порядке:
// М[ row ][ col ] === *( elements + row * stride + col)
// Для удобного доступа к элементам подматриц используется
// вспомогательная структура - дескриптор (под)матрицы: typedef struct {
int width; 11 количество столбцов
int height; // количество строк
int stride; 11 размер одной строки в элементах
float* elements; 11 указатель на массив элементов
} Matrix;
// Для работы с матрицами будет запускаться двумерная сетка
// Размер блока потоков по координатам х и у:
#define BLOCK_SIZE 16
// Определение прототипа ядра
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);
// Функция для GPU, возвращающая элемент подматрицы
__device__ float GetElement( const Matrix A, int row, int col) {
return A.elements[ row * A.stride + col ];
62
// функция записи значения в заданный элемент матрицы
__device__ void SetElement(Matrix A, int row, int col, float value) {
A.elements[ row * A.stride + col ] = value;
}
// Функция формирования дескриптора подматрицы размером
// BLOCK_SIZE * BLOCK_SIZE, которая размещена на col
// столбцов правее и на row строк ниже левого верхнего угла // данной матрицы
__device__ Matrix GetSubMatrix(Matrix A, int row, int col)
{
Matrix Asub;	//	возвращаемый дескриптор
Asub.width	=	BLOCK_SIZE;
Asub.height	=	BLOCK_SIZE;
Asub.stride	=	A.stride;
Asub.elements = &A.elements[ A.stride * BLOCK_SIZE * row + BLOCK_SIZE * col ];
return Asub;
}
// Умножение матриц - код хост-компьютера
// Предполагается, что размеры всех матриц делятся нацело
// на BLOCK_SIZE
void MatMul(const Matrix A, const Matrix B, Matrix R) { // Формирование матриц d_A и d_B на основе аргументов Matrix d_A; //объявление матрицы А d_A.width = d_A.stride = A.width;
d_A.height = A.height;
size_t size = A.width * A.height * sizeof(float);
cudaMalloc(&d_A.elements, size);
cudaMemcpy(d_A.elements, A.elements, size, cudaMemcpyHostToDevice);
Matrix d_B; // объявление матрицы В d_B.width = d_B.stride = B.width;
d_B.height = B.height;
size = B.width * B.height * sizeof(float);
cudaMalloc(&d_B.elements, size);
cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice);
// выделение памяти GPU под матрицу d_R
Matrix d_R; // объявление матрицы R
d_R.width = d_R.stride = R.width;
d_R.height = R.height;
size = R.width * R.height * sizeof(float);
cudaMalloc(&d_R.elements, size);
// Объявление переменных для измерения времени
// работы ядра
float kernelTime;
63
cudaEvent_t beginEvent, endEvent;
cudaEventCreate (&beginEvent);
cudaEventCreate (&endEvent);
cudaEventRecord (beginEvent > 0);
// Запуск ядра:
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
MatMulKernel«< dimGrid, dimBlock »>(d_A, d_B, d_R);
cudaEventRecord (endEvent, 0);
cudaEventSynchronize (endEvent);
// Вычислить разницу моментов времени начала
// и конца работы
cudaEventElapsedTime (&kernelTime, beginEvent, endEvent);
// Так можно вывести на консоль время работы ядра:
// printf(«Kernel time: %.2f milliseconds\n», kernelTime);
// Пересылка матрицы-результата из памяти GPU
// в основную память хост-компьютера
cudaMemcpy(R.elements, d_R.elements, size, cudaMemcpyDeviceToHost);
__ // Обработка полученного произведения матриц:
// Освобождение памяти GPU
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_R.elements);
}
// Собственно программа ядра
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix R)
{
// Координаты строки и столбца блока потоков
int blockRow = blockldx.y;
int blockCol = blockldx.x;
// Каждый блок потоков готовит свою подматрицу результата Matrix Rsub = GetSubMatrix(R, blockRow, blockCol);
// Каждый поток вычисляет один элемент подматрицы путем
// накопления результата в регистровой переменной Rvalue
float Rvalue = 0;
// выборка номеров строки и столбца внутри подматрицы Rsub
int row = threadldx.y;
int col = threadldx.x;
// Цикл по всем подматрицам сомножителей, которые
// требуются для вычисления подматрицы результата
for(int m = 0; m < (A.width / BLOCK_SIZE); ++m) {
// получить очередную подматрицу из А
Matrix Asub = GetSubMatrix(A, blockRow, m);
// получить очередную подматрицу из В
64
Matrix Bsub = GetSubMatrix(B, m, blockCol);
// выделить блоки разделяемой памяти для подматриц _shared_ float As[ BLOCK_SIZE ][ BLOCK_SIZE ]; _shared_ float Bs[ BLOCK_SIZE ][ BLOCK_SIZE ];
// загрузить подматрицы из медленной общей памяти GPU // в быструю разделяемую.
// Каждый поток перепишет ровно по одному элементу
As[ row ][ col ] = GetElement(Asub, row, col);
Bs[ row ][ col ] = GetElement(Bsub, row, col);
// Синхронизировать все потоки внутри одного блока
// для того, чтобы быть уверенным в наличии
// всех элементов подматриц
// в разделяемой памяти перед началом вычислений
__syncthreads();
// Каждый поток вычисляет свой элемент подматрицы
for (int ie =0; ie < BLOCK_SIZE; ++ie)
Rvalue += As[ row ][ ie ] * Bs[ ie ][ col ];
// Синхронизировать потоки перед загрузкой
// двух новых подматриц для следующей итерации
__syncthreads();
} // Конец тела цикла перебора подматриц сомножителей
// Запись вычисленной подматрицы в общую память,
// каждый поток пишет свой вычисленный элемент в результат SetElement(Rsub, row, col, Rvalue);
} // Конец тела ядра
Необходимо отметить, что для краткости из этой программы были удалены все проверки возникновения ошибок при вызове функций библиотеки CUDA-runtime. Каждая вызываемая функция формирует и возвращает код своего завершения (значение типа cudaError_t), который можно получить, вызывая функцию
cudaError_t cudaGetLastError();
Если значение кода завершения не равно CUDA_SUCCESS, то выполнять последующие обращения к функциям библиотеки, как правило, смысла не имеет, поэтому обычно после вызова каждой CUDA-функции пишут фрагмент наподобие такого (переменная err должна быть ранее объявлена как cudaError_t):
if((err = cudaGetLastError()) != CUDA_SUCCESS) { printf(«\nError r%s}} code = %d», err,
cudaGetErrorString(err));
exit(l);
}
В реальные программы эти проверки включать необходимо.
65
4.6.	Библиотека CUDA-runtime, основные функции
В последовательной части этой программы были использованы вызовы некоторых функций библиотеки CUDA-runtime для управления глобальной памятью GPU.
Всего таких функций в библиотеке очень много, для детального изучения следует обращаться к оригинальной документации по CUDA, которая обновляется с каждой новой версией. Здесь будут описаны некоторые часто используемые функции.
4.6.1.	Перед началом работы с графическим процессором необходимо получить его параметры и характеристики, используя вызов функции:
cudaError_t cudaGetDeviceProperties(cudaDeviceProp* props, int device);
Эта функция заполняет структуру, адрес которой содержит указатель props. Используя поля этой структуры, можно узнать, например, объемы глобальной, разделяемой и регистровой памяти (две последние — на один потоковый мультипроцессор), предельные размеры сетки в блоках и блоков в потоках по каждому измерению и другие важные для исполнения программы параметры устройства.
4.6.2.	Для управления глобальной памятью GPU как минимум нужны функции cudaMalloc, cudaFree и cudaMemcpy.
cudaError_t cudaMalloc (void **devPtr, size_t size);
выделяет блок глобальной памяти устройства размером size и заносит адрес этого блока в указатель devPtr. Возвращает либо код успешного завершения CUDA_SUCCESS, либо код ошибки, если блок памяти распределить не удалось.
cudaError_t cudaFree (void *devPtr);
возвращает ранее выделенный блок глобальной памяти в пул свободной памяти. Возвращает либо код успешного завершения, либо код ошибки (например, если этот блок ранее не выделялся).
cudaError_t cudaMemcpy(void *dst, const void *src,
size_t count,
enum cudaMemcpyKind kind);
66
копирует блок памяти размером count байт, расположенный по адресу src, в память по адресу dst. Один из этих блоков размещается в основной памяти компьютера, второй — в глобальной памяти графического процессора. Где именно находятся эти блоки, определяется значением аргумента kind (допустимы значения cudaMemcpyHostToDevice и cudaMemcpyDeviceToHost).
4.6.3.	Для профилирования и отладки CUDA-программ можно использовать функции управления событиями. Все такие функции оперируют со структурами данных типа cudaEvent_t. Как приведено в примере в п. 4.5, вначале структуру-событие следует создать, вызывая функцию
cudaError_t cudaEventCreate (cudaEvent_t* event)
После этого в момент возникновения самого события (например, при запуске ядра на выполнение), нужно вызвать функцию
cudaError_t cudaEventRecord (cudaEvent_t* event
[, cudaStream_t stream ]);
Она копирует (захватывает) в указанное событие состояние указанной вторым аргументом очереди работ (если аргумент опущен, то очереди по умолчанию). Никакие последующие операции с очередью не приведут к изменению сохраненного состояния. Обеспечить завершение всех работ, находящихся в очереди в момент фиксации указанного события можно с помощью функции:
cudaEventSynchronize (cudaEvent_t* event);
Следующая функция позволяет вычислить разницу между моментами фиксации указанных событий (в миллисекундах) и записать ее в переменную, переданную как первый аргумент:
cudaEventElapsedTime ( float* time, cudaEvent_t firstEvent, cudaEvent_t secondEvent);
Множество других полезных функций CUDA-runtime, функций CUDA-driver’a и встроенных функций, выполняемых на графическом процессоре, рекомендуется изучать по документации к используемой версии CUDA [18].
67
Контрольные вопросы и задания
1.	Что такое ядро в терминологии CUDA?
2.	Как задается размерность и количество потоков ядра, выполняемого графическим процессором?
3.	Перечислите виды и характеристики памяти, доступной из программы GPU.
4.	Как указать требуемое размещение переменной в памяти GPU (регистровой, разделяемой, глобальной, памяти текстур,...)?
5.	Перечислите виды и характеристики памяти, доступной из программы CPU при использовании архитектуры CUDA.
6.	Перечислите встроенные векторные типы данных расширения языка С и объясните смысл их наименований.
7.	Как в программе для CPU обеспечить синхронизацию с программой для GPU?
8.	Если к одному CPU подключено несколько видеокарт NVIDIA архитектуры CUDA, то каким образом можно обеспечить их одновременную загрузку?
9.	Перечислите виды и характеристики памяти GPU, доступной из программы основного процессора.
10.	Какие ограничения наложены на функции, которые должны выполняться на GPU?
11.	Что такое мультипроцессор в терминологии CUDA?
12.	Какие ограничения наложены на функции, выполняемые в графическом процессоре?
13.	Как определить версию и технические характеристики графического процессора NVIDIA?
14.	Перечислите и охарактеризуйте группы библиотечных функций Run-time библиотеки CUDA.
15.	Как измеряется время выполнения ядра?
16.	Что такое сетка, блок, поток в терминологии CUDA?
17.	Что такое тип данных dim3?
18.	Из каких компонент состоит потоковый мультипроцессор?
19.	Какова максимальная размерность сетки в блоках?
20.	Перечислите и охарактеризуйте функции компонентов программного обеспечения CUDA.
21.	Какую модель параллелизма реализует архитектура CUDA?
22.	Что делает препроцессор nvcc?
23.	Можно ли получить указатель на функцию, выполняемую графическим процессором?
24.	Что такое локальная память потока?
25.	С помощью какого расширения синтаксиса запускается ядро CUDA?
26.	Перечислите и охарактеризуйте дополнительные виды функций, существующие в расширении синтаксиса языка С для CUDA.
27.	Какие виды памяти доступны из программ для GPU?
68
28.	Как синхронизируются процессы в CPU и в GPU?
29.	Можно ли одновременно использовать технологии ОрепМР и CUDA?
30.	Что такое частная память потока? Каков ее объем?
31.	Перечислите имена полей, используемые во встроенных векторных типах данных.
Тема 5
OPENCL: ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ НА ГЕТЕРОГЕННОМ КОМПЬЮТЕРЕ
OpenCL (Open Computing Language) — это спецификация, описывающая библиотечное расширение языков C/C++, обеспечивающее поддержку параллельных вычислений на совокупностях неоднородного оборудования — гетерогенных компьютерах [19]. Начальная версия OpenCL была разработана компанией Apple, впоследствии для развития и стандартизации спецификации была образована группа разработчиков Khronos Compute [20—22], в которую вошли Apple, NVIDIA, AMD, IBM, Intel, ARM, Motorola и др. Первая версия стандарта OpenCL была опубликована в конце 2008 года.
5.1.	Модель платформы
Платформа OpenCL представляет собой хост-компьютер, с подключенной к нему совокупностью возможно разнотипных устройств [20], для каждого из которых имеется соответствующий программный комплект разработки (Software Development Kit, SDK). Модель платформы показана на рис. 5.1.
SDK существуют для графических процессоров компаний NVIDIA и AMD, для многоядерных процессоров компаний Intel и AMD, для вычислительных акселераторов компании Intel и ряда других. На одном хост-компьютере могут быть установлены и совместно функционировать несколько различных SDK.
Каждое OpenCL-устройство (Compute Device, CD) состоит из одного или нескольких вычислительных блоков (Compute Unit, CU), которые в свою очередь состоят из одного или более элементов-обработчиков (Processing Elements, РЕ).
Устройства в целом и их вычислительные блоки обычно реализуют вычисления согласно принципу ОКМД, но могут работать и согласно принципу ОПМД.
70
5.2.	Модель памяти OpenCL.
В модели памяти стандарта OpenCL [20, 21] предусмотрены четыре различных типа памяти, как показано на рис. 5.2.
Глобальная память (Global Memory). Эта память предоставляет доступ на чтение и запись и хост-компьютеру, и всем элементам-обработчикам любого вычислительного блока OpenCL-устройства. Операции записи и чтения глобальной памяти могут кэшироваться, если соответствующая аппаратура встроена в CD, но кэши «прозрачны» для программиста-разработчика. Программы, выполняемые обработчиками, не могут явно адресовать элементы памяти кэша. Объем глобальной памяти может составлять единицы и десятки гигабайт.
Константная память (Constant Memory). Область глобальной памяти, в которую записывать данные может только хост-компьютер. Обработчики могут только читать данные из этой памяти. Память констант также может кэшироваться. Объем константной памяти обычно значительно меньше, чем глобальной (порядка десятков килобайт).
Локальная память (Local memory). Область памяти, доступная только для элементов-обработчиков данного вычислительного блока как на чтение, так и на запись. Она может ис
71
пользоваться для хранения данных, общих для всех программ, исполняемых обработчиками. Локальная память физически обычно реализуется на чипе OpenCL-устройства и по скорости доступа сравнима с кэшем. Ее объем относительно невелик и составляет десятки килобайт.
Compute Device
Compute Device Memory
Host
Put. 5.2. Модель памяти OpenCL
Частная память (Private memory). Область памяти, принадлежащая исключительно одному элементу-обработчику. Данные, хранящиеся в частной памяти одного обработчика, не видны никаким другим. Частная память, так же как локальная, обычно реализована на чипе и имеет такое же быстродействие. Альтернативно (или частично) она может проецироваться на выделенную область в глобальной памяти устройства. Суммарный объем частной памяти одного вычислительного блока обычно также составляет десятки килобайт.
Спецификация OpenCL определяет 4 уровня памяти, но не накладывает никаких требований на ее реализацию в аппаратуре. Все 4 уровня могут находиться в глобальной памяти, и разделение типов может осуществляться на уровне драйвера.
72
И напротив, может быть реализовано и физическое разделение уровней памяти, обусловленное структурой чипа. Тем не менее, при разработке OpenCL программ следует использовать данные о соотношении объемов и скоростей различных уровней памяти для наиболее эффективного размещения данных.
5.3.	Модель исполнения OpenCL-программ
OpenCL-приложение запускается на хост-компьютере, определяет перечень доступных платформ и устройств, затем формирует и отправляет устройствам запросы на выполнение вычислительных работ, содержащие коды программ и данные для них.
Выполняемая OpenCL-программа текстуально состоит как минимум из двух частей: основная, исполняемая на хост-компьютере, и так называемое ядро (kernel) — специальным образом оформленная функция, предназначенная для исполнения на OpenCL-устройстве [18—21]. Ядер может быть определено несколько. Каждое ядро описывает исполняемый одним элементом-обработчиком вычислительного блока процесс обработки аргументов функции, называемый потоком. Как правило, телом этой функции является тело некоторого цикла последовательного прототипа параллельной программы, т. е. выполнение одной итерации цикла соответствует одному потоку. Ядру должны быть переданы аргументы, определяющие данные для обработки и место в памяти для сохранения результатов выполнения. Все потоки получают одинаковые значения аргументов ядра. Ядром называется функция, «вызываемая» из основной программы хост-компьютера, но выполняющаяся во многих экземплярах на OpenCL-устройстве. Одна функция не может явно возвращать множество результатов своей работы, поэтому тип возвращаемого ею значения всегда должен быть void.
При запуске ядра на выполнение основная часть программы передает системе поддержки вычислений OpenCL параметры так называемого пространства индексов совокупности потоков, называемого NDRange. Пространство индексов может быть 1-, 2- или 3-мерным. Выбор размерности NDRange определяется удобством для конкретного алгоритма: в случае работы с трехмерными моделями и массивами удобно индексировать потоки по трехмерным координатам, в случае работы с изображениями или двумерными сетками — удобнее, когда размерность пространства индексов равна 2, если же решаемая
73
задача имеет линейный характер, то достаточно одномерного NDRange. Каждый поток имеет свои собственные координаты в этом пространстве и может получить их, используя встроенные функции OpenCL.
OpenCL-программа может быть запущена на хост-компьютере, в состав которого входит графический процессор. Поэтому при написании функции-ядра необходимо учитывать возможность разбиения NDRange на блоки в соответствии с физической организацией OpenCL-устройства. Согласно терминологии OpenCL блок потоков называется Work-group. Work-group выполняется одним CU. Каждый CU, в свою очередь, состоит из нескольких обработчиков. Вычислительный поток, выполняемый одним обработчиком, называется Workitem.
Каждый CU работает как ОКМД-система и имеет свою собственную память наряду с основной памятью устройства. Группа потоков, исполняемая одним CU, имеет свои групповые индексы в NDRange, а каждый ее поток — локальные индексы внутри группы (и может их получить наряду с глобальными индексами). Пример возможной организации пространства индексов NDRange приведен на рис. 5.3.
Gx = 15	-----Т~
Wx = 3	Lx~ 5
Рис 5.3. Один из возможных вариантов пространства индексов
74
На этом примере показано двумерное пространство групп потоков размером 3x3. Каждая группа имеет размерность 5 на 5 потоков, поэтому общие размерности пространства индексов по координатам х и у равны 15, всего — 225 потоков. Каждый поток имеет свои индексы и в группе (локальные), и в пространстве NDRange (глобальные). При запуске этой совокупности потоков она вся будет выполняться одним OpenCL-устройством. Разные группы этой совокупности, возможно, будут выполняться разными вычислительными блоками, но каждая группа — всегда одним, и при этом одним и тем же, блоком. Другими словами, группы потоков не могут мигрировать с одного вычислительного блока на другой, т. е. не могут начинаться на одном, продолжаться на другом и заканчиваться на третьем.
5.4.	Модель структуры OpenCL-программ
В OpenCL (аналогично CUDA) программа-приложение разделяется на две части: первая часть — управляющая (последовательная), вторая — вычислительная (параллельная). Однако, в отличие от CUDA, в которой компиляция и той, и другой частей выполняется единовременно, в OpenCL управляющая часть компилируется заблаговременно, а вычислительная — в процессе исполнения управляющей частью после выбора того устройства, на котором она будет запускатья.
Это обусловлено ориентацией стандарта OpenCL на поддержку различных, в том числе — возможно еще не существовавших на момент написания приложения устройств. Ясно, что компиляция и создание исполняемого кода для еще несуществующего устройства невозможна. Поэтому текст функции-ядра, т. е. вычислительной части программы, включается в управляющую часть в «исходном» виде, т. е. как символьная строка. Возможен и достаточно часто используется вариант, при котором текст функции-ядра оформляется в виде отдельного файла, который считывается управляющей частью непосредственно перед компиляцией средствами выбранного SDK. Исходный текст функции-ядра компилируется средствами OpenCL непосредственно в процессе работы приложения для выбранного в данный момент вычислительного устройства. Это выполняется при каждом запуске OpenCL-программы.
75
Тем не менее в OpenCL предусмотрена возможность сохранить скомпилированный код ядра в бинарном представлении и загружать его впоследствии каждый раз при инициализации вычислительной программы. Это позволит сэкономить некоторое время (загрузка готового бинарного кода выполняется быстрее загрузки компилятора и его работы), однако приведет к потере универсальности, поскольку этот код будет работать только на одном типе устройств.
Общая структура управляющей части OpenCL-программы выглядит более сложной, чем структура аналогичной CUDA-программы. В OpenCL необходимо выполнять дополнительные действия по конфигурированию среды выполнения и подготовке кода для устройства. Среда выполнения в терминах OpenCL именуется контекстом, это структура данных, связанная с платформой (платформа — это совокупность однотипных устройств одного производителя, по сути — другое название для SDK), вычислительными устройствами, исходными и бинарными кодами программ параллельной части, очередями работ к устройствам, блоками памяти для хранения обрабатываемых данных.
Для создания контекста и управления параллельной частью приложения последовательная часть вызывает функции библиотеки OpenCL. В имени каждой функции есть общий префикс cl. Каждая OpenCL-функция формирует код завершения, только значение CL_SUCCESS свидетельствует о возможности продолжения работы. Любое другое значение — это код возникшей ошибки. Код завершения необходимо проверять после каждого вызова любой функции (как и при написании CUDA-программ) и организовывать реакцию на возможные ошибки.
Типичная последовательность этапов выполнения OpenCL-приложения может выглядеть так (для каждого этапа в скобках приводятся имена OpenCL-функций, используемых при его выполнении).
1.	Получить общую и детальную информацию о платформах и устройствах (clGetPlatformIDs(), clGetPlatformlnfoQ), clGet DevicelDsO, clGetDevicelnfofX).
2.	Выбрать нужное устройство (устройства) и создать для него (для них) контекст (.clCreateContextO').
3.	Создать очередь работ устройства (clCreateCommand-Queue (), clCreateCommandQueueWithProperties ()).
4.	Создать код ядра для выбранного устройства из текста функции iclCreateProgramWithSourceO, clBuildProgram(J, clCreateKemel () ).
76
5.	Выделить память для данных на устройстве (clCreateBujferQ').
6.	Скопировать данные из основной памяти хост-компьютера в глобальную память устройства (clEnqueueWriteBujferQ').
6.	Привязать параметры к функции-ядру (clSetKernelArgO).
7.	Поставить ядро в очередь к устройству (clEnqueueND RangeKernelXX).
8.	Скопировать результаты работы ядра из глобальной памяти устройства в основную память хост-компьютера (clEnqueueReadBufferO').
9.	Освободить занятые ресурсы (clReleaseMemObjectO, clRelease KernelO, clReleaseProgramQ, clReleaseContextO, clReleaseCommand Queue O').
10.	Обработать результаты выполнения ядра
11.	Возможно — повторить шаги А—10 (или 5—10, 6—10, 7—10).
12.	Завершить работу приложения.
5.5. Пример OpenCL-программы
В качестве примера разработки приложения на OpenCL приведем программу умножения матриц, до некоторой степени аналогичную примеру из темы 4. Чтобы эта программа могла работать, нужно, чтобы в файле с именем «MatrixMulKernel.cl» было записано ядро вида:
// максимальный размер группы потоков по координатам х и у:
#define WGS 32
___kernel void mMul(const __global int* mSizes, // wA,hA, wB, hB
const __global float* A,
const __global float* B, _global float* R)
{
__local float subA[ WGS ][ WGS ];
// для фрагмента матрицы A
__local float subB[ WGS ][ WGS ];
// для фрагмента матрицы В
// глобальные идентификаторы потока:
const gid0 = get_global_id(0).; // по измерению 0
const gidl = get_global_id(l).; // по измерению 1
// фактический размер группы потоков
// (по каждому измерению):
const int wgSize = get_local_size( 0 );
if( wgSize == get_local_size( 1 )) {
// если равны - работаем
77
// индексы потока в группе
const int row = get_local_id( 0 );
const int col = get_local_id( 1 );
// эти размеры матриц используются неоднократно:
const int wA = mSizes[0];
const int wB = mSizes[2];
// Индексы фрагмента матрицы-результата,
// вычисляемого группой потоков
const int gr0 = get group id( 0 ) * wgSize;
const int grl = get_group_id( 1 ) * wgSize;
float value = 0.0f; // накопительная переменная потока
// Перебор фрагментов множимого и множителя:
for ( int iF=0; iF < wA; iF += wgSize ) {
// Загрузка очередных фрагментов в разделяемую память subA[ row ][ col ] = *(А + (gr0 + row)* wA + iF + col); subB[ row ][ col ] = *(B + (iF + row)* wB + grl + col); // Ожидание завершения загрузки всеми потоками группы barrier( CLK_LOCAL_MEM_FENCE );
// Пополнение накопительной переменной:
for ( int k=0; k < wgSize; k++ ) {
// проверка нахождения в пределах фрагмента:
if((k + iF) < wA)
value += subA[ row ][ к ] * subB[ к ][ col ];
}
// Ожидание завершения использования
//текущих фрагментов
barrier( CLK_LOCAL_MEM_FENCE );
}
// Проверка нахождения в пределах матрицы-результата: if((gid0 < wB) && (gidl < mSizes[l])) {
R[ gid0 * wB + gidl ] = value; // сохраняем результат }
}
// Если размеры группы по измерениям 0 и 1 различаются:
else if((gid0 == 0) && (gidl == 0)) // только один поток printf(«\nKernel error: group width (=%d) ang group» +
" height (=%d) must be the same.\n No result computing.",
wgSize, get_local_size(l));
}
Последовательная часть OpenCL-программы, запускающей это ядро, выглядит так:
int main(int argc, char** argv[]) {
int mSizes[] = { 324, 219, 219, 324 };
int sizeA = mSizes[0] * mSizes[l];
int sizeB = mSizes[2] * mSizes[3];
78
int sizeR = mSizes[l] * mSizes[2];
char *mMult;
FILE *fp;
size_t program_size;
fp = fopen("MatrixMulKernel.cl", "r");
fseek(fp, 0, SEEK_END); program_size = ftell(fp) + 100; rewind(fp); // или fseek( fp, 0, SEEK_START );
mMult = (char*)malloc(program_size + 1);
memset(mMult, 0, program_size);
// Считывание ядра OpenCL из файла kernel.cl fread(mMult, sizeof(char), program_size, fp); fclose(fp);
// Определение необходимых переменных и массивов: size_t sizeWork[2];
size_t sizeWG[2];
float *h_a, *h_b, *h_r;
// Указатели на массивы в памяти CPU cl_int clem;
// Для получения кода возврата OpenCL функций
cl_uint qty_platforms = 0; // Количество OpenCL-платформ // Список идентификаторов OpenCL-платформ: cl_platform_id* platforms;
// Количество OpenCL-устройств для каждой платформы: cl_uint *qty_devices;
// Массив идентификаторов OpenCL-устройств: cl_device_id **devices;
cl_context context; // Дескриптор контекста cl_command_queue qE; // Дескриптор очереди команд cl_program program; // Дескриптор программы-ядра cl_kernel kernel; // Дескриптор объекта-ядра // Указатели на данные в памяти OpenCl-устройства, // они же аргументы ядра cl_mem arg0;
cl_mem argl;
cl_mem arg2;
cl_mem arg3;
// Дескриптор события, которое будет связано с ядром: cl_event kEvent;
// Переменные для отсчетов времени: long long tStart = 01;
long long tEnd = 01;
// 1. Получение количества доступных платформ.
clem = clGetPlatformIDs(0, NULL, &qty_platforms);
if (qty_platforms == 0) {
fprintf(stderr, "No Open-CL-devices\n"); exit(-l);
}
// Обработка ошибок. Далее аналогичные блоки опущены
79
if (clerr != CL_SUCCESS) {
fprintf(stderr, "Error, code = %d.\n”, clerr);
exit(-l);
}
// Выделение памяти для хранения информации о конфигурации
// и параметрах устройств
platforms = (cl_platform_id*)malloc(
sizeof(cl_platform_id) * qty_platforms);
devices = (cl_device_id**)malloc(
sizeof(cl_device_id*) * qty_platforms);
qty_devices = (cl_uint*)malloc(sizeof(cl_uint) * qty_ platforms);
// Получение списка идентификаторов платформ
clerr = clGetPlatformIDs(qty_platforms, platforms, NULL);
for (ui = 0; ui < qty_platforms; ui++) {
// Получение количества имеющихся OpenCL-устройств
// для каждой платформы
clerr = clGetDeviceIDs(platforms[ui],
CL_DEVICE_TYPE_ALL, 0, NULL, &qty_devices[ui]);
// Получение списка идентификаторов OpenCL-устройств
if (qty_devices[ui]) {
devicesfui] = (cl_device_id*)malloc( qty_devices[ui] * sizeof(cl_device_id));
clerr = clGetDeviceIDs(platforms[ui], CL_DEVICE_TYPE_ ALL, qty_devices[ui], devices[ui], NULL);
}
}
// Получение максимального размера группы потоков:
unsigned char retVal[9];
size_t actualsize;
clGetDeviceInfo(devices[0][0], CL_DEVICE_MAX_WORK_GROUP_ SIZE,
8, retVal, &actualSize);
int val = *(int*)retVal;
// Установка размеров группы потоков:
switch (val) {
case 256: case 512:
sizeWG[0] = 16; sizeWG[l] = 16; break;
case 1024: case 2048:
sizeWG[0] = 32; sizeWGfl] = 32; break;
default:
printf("\nCanJt use this device")^ exit(-l);
}
// Установка размеров сетки потоков:
sizeWork[0] = mSizes[l];
if ((mSizes[2] % sizeWG[0]) != 0)
sizeWork[0] += sizeWG[0] - mSizes[2] % sizeWG[0];
80
sizeWorkfl] = mSizes[2];
if ((mSizesfl] % sizeWG[l]) != 0)
sizeWork[l] += sizeWG[l] - mSizes[l] % sizeWG[l];
// 2. Формирование контекста. Создается контекст,
// только для платформы 0 (это упрощение,
//реально нужно для всех).
clem = CL_SUCCESS;
context = clCneateContext(0, qty_devices[0],
devices[0], NULL, NULL, &clerr);
// 3. Создание очереди команд для устройства 0 платформы 0.
qE = clCreateCommandQueue(
context, devices[0][0], CL_QUEUE_PROFILING_ENABLE, &clerr);
// 4. Компиляция программы для OpenCL-устройства.
program = clCreateProgramWithSource(
context, 1, (const char**)&mMult, NULL, &clerr);
clem = clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
// Обработка возможных ошибок компиляции - это важно!
if (clem = CL_BUILD_PROGRAM_FAILURE) {
// получить размер лог-файла ошибок компилятора: size_t log size;
clGetProgramBuildInfo(program, devices[0][0], CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
char *log = (char *)malloc(log_size);
// получить лог-файл сообщений об ошибках в ядре
clGetProgramBuildInfo(program, devices[0][0], CL_PROGRAM_BUILD_LOG, log_size, log, NULL);
// вывести сообщения об ошибках для пользователя fprintf(stderr, "%s\n", log);
exit(-l); // завершить работу приложения
}
// Декларация ядра в программе
kernel = clCreateKernel(program, "mMul", &clerr);
// Резервирование памяти в памяти хост-машины
h_a = (float *)malloc(sizeA * sizeof(float));
h_b = (float *)malloc(sizeB * sizeof(float));
h_r = (float *)malloc(sizeR * sizeof(float));
memset(h_r, 0, sizeR * sizeof(float));
// 10. Формирование значений умножаемых матриц.
for (int i = 0; i < sizeA; i++) {
h_a[i] = (i + mSizes[0] - 1) / 32;
}
for (int i = 0; i < sizeB; i++) { h_b[i] = i % (mSizes[2] + 1);
}
// 5. Формирование аргументов для ядра
arg0 = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR, sizeof(mSizes) * sizeof(size_t), mSizes, NULL);
81
argl = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR, sizeA * sizeof(float), h_a, NULL); arg2 = clCreateBuffer(context, CL_MEM_READ_ONLY |
CL_MEM_COPY_HOST_PTR, sizeB * sizeof(float), h_b, NULL); arg3 = clCreateBuffer(context, CL_MEM_WRITE_ONLY |
CL_MEM_COPY_HOST_PTR, sizeR * sizeof(float), h_r, NULL); // 6. Передача аргументов в ядро. clSetKernelArg(kernel, 0, sizeof(cl_mem), (void*)&arg0); clSetKernelArg(kernel, 1, sizeof(cl_mem), (void*)&argl); clSetKernelArg(kernel, 2, sizeof(cl_mem), (void*)&arg2); clSetKernelArg(kernel, 3, sizeof(cl_mem), (void*)&arg3); // Инициализация дескриптора события kEvent = clCreateUserEvent(context, &clerr);
// 7. Запуск ядра на вычисления.
clerr = clEnqueueNDRangeKernel(qE, // queuel - очередь kernel,	// Запускаемое ядро
2,	// Количество измерений сетки потоков
NULL, // В стандарте пока не задействовано sizeWork, // Размерности сетки по ее измерениям sizeWG, // Размерности группы потоков по ее измерениям 0, NULL, &kEvent);
if (clerr != CL_SUCCESS) {
fprintf(stderr, "\nOpenCl error, code: %d", clerr, ); exit(-l);
}
// 8. Копирование результатов в основную память хостмашины. clerr = clEnqueueReadBuffer(qE, arg3, CL_TRUE, 0, sizeR * sizeof(float), h_r, 0, NULL, NULL);
// Получение данных о моментах начала и конца работы ядра clerr = clGetEventProfilingInfo(
kEvent, CL_PROFILING_COMMAND_START,
sizeof(tStart), &tStart, NULL);
clerr = clGetEventProfilingInfo(
kEvent, CL_PROFILING_COMMAND_END, sizeof(tEnd), &tEnd, NULL);
// Вывод на консоль времени работы ядра printf("\nKernel time:%lld”, tEnd - tStart); // 9. Освобождение памяти в OpenCL-устройстве и ресурсов. clReleaseMemObject(argl);
clReleaseMemObject(arg2); clReleaseMemObject(arg3); clReleaseKernel(kernel); clReleaseProgram(program); clReleaseCommandQueue(qE); clReleaseContext(context); free(qty_devices);
free(devices); free(platforms);
82
// 10. Сравнение результата с эталоном.
int errCnt = 0;
int colsA = mSizes[0], rowsA = mSizesfl],
colsB = mSizes[2], rowsB = mSizes[3].;
for (int i = 0; i < rowsA; i += 1) {
for (int j = 0; j < colsB; j += 1) {
float curE = 0.0f;
for (int к = 0; к < colsA; к += 1)
curE += *(h_a + i * colsA + k) * *(h_b +
к * colsB + j);
float e = curE - *(h_r + i * rowsA + j);
if ((e > 0.001f) || (e < -0.001f)) {
if (errCnt < 20)
printf("\ni=%d,j=%d, R=%f, e=%f", i, j,
*(h_r + i * rowsA + j), curE);
errCnt += 1;
}
}
}
// 11. Продолжение работы, возможно,
// повторение запуска ядра
// 12. Завершение работы return 0;
}
Необходимо отметить, что приведенная программа имеет демонстрационный характер, в частности — в ней используются специально сгенерированные матрицы. Размеры перемножаемых матриц указаны в первой строке программы. Как видно из этой строки, размеры матриц не обязательно должны быть кратными размерам групп потоков.
Множество других полезных функций системы поддержки OpenCL и встроенных функций, выполняемых на Compute Device, рекомендуется изучать по документации к используемой версии OpenCL [22].
Контрольные вопросы и задания
1.	Что такое поток, группа потоков и пространство NDRange?
2.	С помощью какой функции можно получить параметры OpenCL-устройства?
3.	Как можно измерить время исполнения ядра?
4.	Сколько очередей работ может быть создано для одного OpenCL-устройства?
83
5.	Как задается размерность и количество потоков ядра, выполняемого OpenCL-устройством?
6.	С помощью каких функций осуществляется управление памятью OpenCL-устройства из хост-программы?
7.	Когда нужна и когда не нужна синхронизация потоков в группе?
8.	Как указать требуемое размещение переменной в памяти OpenCL-устройства (частной, локальной, глобальной)?
9.	Что такое контекст OpenCL-устройства?
10.	Какую модель вычислений обычно реализует OpenCL-устройство?
11.	Что такое встроенные функции? Для чего они используются?
12.	Если при компиляции ядра были обнаружены ошибки, то как получить сообщения о них?
13.	С помощью каких функций выполняется привязка значений аргументов ядра?
14.	Как в хост-программе следует обрабатывать возможные ошибки управления OpenCL-устройством?
15.	Что такое очередь работ OpenCL-устройства?
16.	Можно ли синхронизировать все потоки группы? А все потоки одного исполнения ядра?
17.	Каким может быть количество измерений пространства NDRange?
18.	Что такое события? Как и для чего можно их использовать?
19.	Что такое локальная память группы потоков? Как ее можно эффективно использовать?
20.	В какой момент выполняется компиляция ядра?
21.	Перечислите основные этапы подготовки к запуску ядра на исполнение.
22.	Что такое OpenCL-платформа? Как узнать количество платформ, доступных на данном хост-компьютере?
23.	Перечислите виды памяти OpenCL-устройства, доступной из программы хост-программы.
24.	Как можно синхронизировать потоки в группе?
25.	Перечислите и охарактеризуйте функции компонентов программного обеспечения OpenCL.
26.	Как поток может узнать свои координаты в группе и в общем пространстве NDRange?
27.	Что такое частная память потока? Как узнать ее объем?
Тема б
MPI: ПАРАЛЛЕЛЬНЫЕ ВЫЧИСЛЕНИЯ НА СИСТЕМАХ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ
Аббревиатура MPI расшифровывается так: Message passing interface (интерфейс передачи сообщений). Этот интерфейс в настоящее время в параллельном программировании для вычислительных систем с распределенной памятью является стандартом де-факто [23, 24], поддерживается Аргоннской национальной лабораторией в США и развивается форумом разработчиков MPIForum [25]. Известны и популярны две основные версии стандарта — 1.1 и 2.0, для которых существует довольно большое количество реализаций интерфейса MPI, как свободно распространяемых, так и коммерческих. В настоящем пособии рассматривается свободно распространяемая реализация MPICH. Буквы СН являются сокращением от слова chameleon, которым разработчики обозначили многоплатфор-менность, т. е. возможность использования на различных платформах — Linux, Windows, MacOS. В MPICH поддерживается вся функциональность версии 1.1 стандарта MPI и некоторые возможности версии 2.0.
Реализация интерфейса MPICH [26] включает в себя следующие составные части:
•	Smpd — демон (в терминах Linux) или служба (в терминах Windows). Постоянно работает на каждом узле сети, на котором возможен запуск/выполнение параллельных MPI-программ. Обслуживает запросы на запуск/останов этих программ и запросы от их ветвей на взаимодействие. Экземпляры MPI-программ, запущенные демоном, называются ветвями;
•	библиотека функций MPI — обеспечивает интерфейс между ветвью параллельной программы и демоном smpd. Для использования функций MPI к программе нужно подклю
85
чать заголовочные файлы (имеющие расширения .h или .hpp) и файлы собственно библиотеки (с расширениями .а, .о, .so для Linux, Jib для Windows);
•	утилиты, с помощью которых осуществляется запуск параллельных программ, настройка среды и собственно MPL Для запуска используется утилита MPIRun (в Windows — MPIexec, имеющая графический интерфейс пользователя — WMPIexec). Для настройки MPI и окружающей среды используются утилиты MPIConfig и MPIRegister;
•	библиотека МРЕ — обеспечивает сбор, обработку и визуализацию результатов обработки статистики взаимодействий параллельных ветвей. Используется для профилирования, отладки и оптимизации параллельных программ.
6.1.	Основные понятия интерфейса MPI
Отличительная особенность программ, написанных с использованием стандарта MPI-1, состоит в том, что в них допускается только статическое распараллеливание, т. е. количество параллельных процессов (ветвей) во время запуска и исполнения программы фиксировано.
Стандарт MPI-1 включает в себя следующие группы функций.
1.	Инициализация, завершение, определение окружения.
2.	Передача сообщений типа «точка-точка».
3.	Коллективные операции взаимодействия.
4.	Создание и использование производных типов данных.
5.	Управление группами ветвей и коммуникаторами.
6.	Управление виртуальными топологиями.
Основными понятиями, на которых базируется вся функциональность интерфейса MPI, являются понятия коммуникатора и группы параллельных ветвей. Любые взаимодействия между ветвями параллельной программы реализуются с обязательным указанием коммуникатора.
Коммуникатором называется указатель (дескриптор) на внутреннюю структуру библиотеки MPI, хранящую сведения о коллективе параллельных ветвей. В MPI понятия коллектива ветвей и группы ветвей различаются. Группа ветвей как объект может существовать самостоятельно, коллектив — это всегда коммуникатор. Имена встроенных коммуникаторов обычно
86
начинаются с префикса МР1_С0ММ_, с него же рекомендуется начинать имена создаваемых в приложении коммуникаторов (интерфейс MPI предоставляет такую возможность).
Коммуникатор коллектива, который включает в себя все ветви приложения, создается автоматически при инициализации библиотеки и называется MPI_COMM_WORLD. В этот момент создаются также:
•	коммуникатор MPI_COMM_SELF. В этот коллектив входит единственная ветвь, вызвавшая функцию инициализации;
•	пустой коммуникатор MPI_COMM_NULL, не содержащий ни одной ветви. Использование этого коммуникатора в вызове любой функции MPI приведет к аварийному завершению ветви, вызвавшей функцию, и аварийному завершению всего приложения.
В программе можно явно создавать сколько угодно новых коммуникаторов, указывая нужные коллективы параллельных ветвей (предварительно нужно создавать соответствующие группы ветвей MPI для указания их в функциях создания коммуникаторов). Любая ветвь может быть членом произвольного количества коммуникаторов.
Нужно помнить, что каждый коммуникатор, как и любой другой объект, создаваемый средствами библиотеки, расходует ресурсы, в частности, оперативную память.
Для приема/передачи сообщений важным является понятие буфера данных. Во многих функциях MPI должен быть указан буфер, содержащий передаваемые данные или предназначенный для принимаемых от других ветвей данных. Каждый буфер всегда определяется тремя параметрами:
1)	адрес (указатель на первый байт буфера);
2)	длина (размер буфера в единицах длины типа, заданного третьим параметром);
3)	тип (в MPI определен набор типов данных, указываемых с помощью символических имен; существует возможность создания собственных, или производных, типов, указываемых с помощью имен переменных, получающих в качестве значений дескрипторы пользовательских типов).
Перечень базовых (предопределенных) типов данных MPI содержится в табл. 6.1.
Важную роль в MPI играет понятие тега. Тегом называется целочисленное значение, используемое для установления соответствия внутри потоков сообщений, передаваемых от одной
87
ветви другой. С помощью тегов решается такая проблема: если ветвь А асинхронно (без синхронизации) отправляет два (или более) сообщения одинакового размера, адресованных ветви В, то в общем случае нет гарантии, что сообщение, отправленное первым, будет первым и доставлено.
Таблица 6.1
Тип данных MPI	Тип данных языка С	Тип данных языка C++	Тип данных языка Fortran
MPI.CHAR	signed char	signed char	Character(l)
MPI.WCHAR	wchar.t	wchar.t	Character(l)
MPI.SHORT	signed short int	signed short int	Integer * 2
MPIJNT	signed int	signed int	Integer * 4
MPI.LONG	signed long int	signed long int	Integer * 4
MPI_UNSIGNED_ CHAR	unsigned char	unsigned char	Character
MPI_UNSIGNED_ SHORT	unsigned short int	unsigned short int	Integer * 2
MPI.UNSIGNED	unsigned int	unsigned int	Integer * 4
MPI_UNSIGNED_ LONG	unsigned long int	unsigned long int	Integer * 4
MPIBOOL (MPI_ LOGICAL)	signed int	bool	Logical
MPI.FLOAT	float	float	Real * 4
MPI.DOUBLE	double	double	Double precision (Real * 8)
MPI_LONG_DOUBLE	long double	long double	Double precision (Real * 8)
MPI.COMPLEX	—	Complex <float>	Complex
MPI_DOUBLE_ COMPLEX	—	Complex <double>	Double complex
MPI_LONG_ DOUBLE. CO MPLEX	—	Complex <long double>	—
88
Окончание табл. 6.1
Тип данных MPI	Тип данных языка С	Тип данных языка C++	Тип данных языка Fortran
MPI_BYTE (8-битный байт)	—	—	—
MPI_PACKED	—	—	—
По самым разным причинам может оказаться, что второе сообщение будет доставлено раньше первого. Если оно будет обработано ветвью В как первое, то результаты работы программы в целом будут неверны. Указание разных значений тегов для этих сообщений как на передающей, так и на приемной стороне позволяет избежать путаницы при доставке сообщений.
Рекомендуется теги объявлять с помощью оператора препроцессора #define, выбирая семантически значимые имена для целочисленных констант, например:
#define FIRST_VECTOR 100
#define COLUMN 101
Для сообщений, обмен которыми выполняется внутри циклических участков программы, в качестве тега рекомендуется использовать индекс соответствующего цикла.
Подавляющее большинство функций MPI возвращают целое значение, равное MPI_SUCCESS, если выполнение функции завершилось без ошибок, и код ошибки в противном случае. Коды ошибок зависят от реализации библиотеки и содержатся в справочных руководствах.
Далее рассматриваются основные группы функций [27—29] и приводится краткое описание функций.
6.2.	Функции инициализации и определения окружения
int MPI_Init(int *argc, char **argv); — до вызова этой функции не может быть вызвана никакая другая функция библиотеки MPI (за исключением функции MPIJ.nitializedQ').
int MPI_FinalizeO; — завершение: после вызова этой функции больше нельзя вызывать какие бы то ни было функции
89
MPI, кроме MPI_initialized (однако в некоторых реализациях допускается повторная инициализация). Для нормального завершения параллельной программы эту функцию должны вызвать все ее ветви.
int MPI_Comm_size(MPI_Comm comm, int *size); — получение количества ветвей в коллективе коммуникатора сотт.
int MPI_Comm_rank(MPI_Comm comm, int *rank); — получение индекса (номера) данной ветви в коммуникаторе сотт.
int MPI_Get_j)rocessor_name (.char *name, int *resultlen); — получение имени узла сети/кластера, на котором выполняется данная ветвь.
double MPI_Wtime(); — получение текущего времени, которое отсчитывается от момента запуска программы.
double MPI_Wtick(); — получение величины тика таймера узла в секундах.
int MPI_Buffer_attach(void ^buffer, int size); — передача библиотеке MPI буфера памяти для некоторых операций обмена сообщениями.
int MPI_Buffer_detach(void *buffer_addr, int size); — извещение библиотеки MPI о том, что указанный буфер памяти использовать больше нельзя. Возврат из функции будет выполнен только тогда, когда все буферизованные операции передачи данных, связанные с этим буфером, будут завершены, т. е. когда этот буфер будет полностью освобожден.
6.3.	Обмен сообщениями типа «точка-точка»
Эта группа функций обеспечивает передачу/прием сообщений между парами ветвей одного коммуникатора. Передача/ прием могут быть блокирующими (выполнение ветви останавливается вплоть до момента завершения обмена данными) и неблокирующими (ветвь получает управление сразу после того, как библиотека обработала аргументы функции).
int MPI_Send(yoid *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm); — передача сообщения с блокировкой. Здесь и далее три аргумента buf, count и datatype полностью определяют буфер с данными. Аргументы dest и tag однозначно определяют индекс ветви-получателя сообщения в коммуникаторе сотт и порядковый номер сообщения.
90
int MPI_Recv(yoid *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status ^status'); — прием сообщения, отправляемого функцией MPI_Send, с блокировкой.
Смысл аргументов функции точно такой же, как и у функции MPIJSend (аргумент source кроме очевидного номера ветви-отправителя в коммуникаторе сотт может иметь значение МР1_ ANY_SOURCE, т. е. отправитель —любая ветвь коммуникатора), однако здесь есть дополнительный аргумент status, являющийся указателем на структуру, поля которой позволяют получить дополнительные сведения о принятом сообщении (их перечень определен исходя из нужд функций без блокировки):
•	int cancelled — содержит 0, если прием завершился успешно, не 0 в противном случае;
•	int count — размер буфера для принимаемых данных (не фактическое количество принятых данных);
•	int MPI_ERROR — код ошибки, если прием завершился с ошибкой (может не совпадать с тем, что возвращает функция MPI_Recv);
•	int MPI_SOURCE — идентификатор ветви-отправителя данных;
•	int MPI_TAG — тег принятого сообщения.
Если ни одна ветвь не отправит в данную ветвь сообщения с указанным тегом, то ветвь-получатель никогда не получит управления, т. е. параллельная программа просто зависнет.
Операции MPI_Send и MPI_Recv являются основными для обменов типа «ветвь-ветвь». Однако их использование в реальных программах сопряжено с известными трудностями и в большинстве случаев приводит к невысокой производительности.
Поэтому в MPI дополнительно поддерживаются следующие разновидности операций передачи.
•	Стандартная передача. Считается завершенной, как только сообщение отправлено по сети, независимо от того, получено оно или нет. Передача сообщения может начинаться, даже если принимающая ветвь еще не вызвала функцию приема.
•	Синхронная передача. Не завершается до тех пор, пока не будет завершен прием сообщения.
•	Буферизированная передача. Вызов функции завершается сразу же, потому что сообщение при этом копируется в системный буфер и там ожидает пересылки.
•	Передача по готовности. Собственно передача начинается только в том случае, если начат прием сообщения. В этот
91
момент без ожидания окончания приема выполняется возврат в ветвь.
Кроме того, каждая из этих четырех разновидностей передачи может быть выполнена в блокирующей или неблокирующей форме (блокирующий прием/передача приостанавливает ветвь на время приема/передачи сообщения). Таким образом, всего существует восемь разновидностей операций передачи, включая уже описанную функцию MPI_Send. В MPI принято соглашение об именах процедур, позволяющее легко определять тип используемой операции по имени функции, которое строится по следующему правилу:
MPI_[I][R|S|B]Send
где I (Immediate) — обозначает неблокирующую операцию; R (Ready) — передача по готовности; S (Synchronous) — синхронный; В (Buffer) — буферизированный (только первая буква после подчеркивания является прописной, все остальные — строчные, как, например, MPI_Ibsend).
Аргументы всех блокирующих (без первой после символа подчеркивания буквы Г) функций точно такие же, как у функции MPI_Send. Добавочным аргументом неблокирующих функций MPI_I*send является дескриптор так называемой квитанции, который в программе должен быть объявлен так: MPI_Request request; (имя переменной может отличаться от имени request). Квитанции используются библиотекой MPI для синхронизации работы ветви и процессов доставки сообщения.
Кроме блокирующей функции приема сообщения MPI_Recv существует ее неблокирующий вариант MPI_Irecv с добавочным аргументом-квитанцией.
Запущенную ранее неблокирующую операцию приема/ передачи в любой момент до ее завершения можно отменить, вызвав функцию
int MPI_Cancel(MPI_Request *request).;
Статус принимаемого сообщения можно получить, не читая из MPI самого сообщения. Это делается либо с помощью блокирующей функции
int MPI_Probe(int source, int tag, MPI_Comm comm, MPI_Status *status);
92
либо с помощью неблокирующей функции
int MPI_Iprobe(int source, int tag, MPI_Comm comm, int flag, MPI_Status *status);
Проверить, завершилась ли одиночная неблокирующая операция, можно путем вызова функции
int MPI_Test(MPI_Request *request, int flag, MPI_ Status *status).;
Эта функция проверяет состояние операции, связанной с квитанцией request, и формирует значение флага false, если операция еще не завершилась, и значение true, если операция завершилась (в этот момент становятся актуальными поля структуры status').
Дождаться завершения одной операции можно, вызвав функцию
int MPI_Wait (MPI_Request *request, MPI_Status *status);
Вместо адреса структуры status может быть указано значение MPI_STATUS_IGNORE. Тогда определить, нормально ли завершилась операция, можно только по коду завершения, возвращаемому функцией MPI_Test или MPI_Wait. Эта функция возвращает управление в ветвь только тогда, когда завершается неблокирующая операция MPI, в вызове которой был указан данный request. В этот момент формируются/модифицируются поля структуры status (если она указана).
Если было запущено более одной неблокирующей операции приема/передачи, то с помощью функции:
int MPI_Testany(int count, MPI_Request array_of_ requests[], int *index, int flag, MPI_Status *status).;
можно узнать, не завершилась ли хотя бы одна из них (любая, ее порядковый номер будет записан в аргумент index, если аргумент flag будет после вызова иметь значение true).
Проверить завершение произвольного количества (но не обязательно всех) операций приема/передачи, связанных с квитанциями из массива array_of_requests, можно с помощью функции
int MPI_Testsome(int incount,
MPI_Request array_of_ requests[],
93
int *outcount, int array_of_indices[], MPI_ Status array_of_statuses[]);
Аргумент incount (как и count в предыдущей функции) определяет размер массивов array_of_requests, array_of_indices и array_of_ statuses, аргумент outcount предназначен для возврата количества завершившихся операций, индексы соотве-ствующих квитанций будут записаны в массив array_of_indices.
Проверить факт завершения всех операций приема/переда-чи, квитанции которых содержатся во втором аргументе, можно с помощью функции
int MPI_Testall(int count, MPI_Request array_of_ requestsf], int flag, MPI_Status array_of_statuses[]).;
Проверить, отменена или нет операция приема/передачи, можно с помощью функции, устанавливающей flag:
int MPI_Test_cancelled(MPI_Status *status, int flag);
Ожидать завершения любой ранее запущенной операции передачи/приема, квитанции на которые указаны в массиве, передаваемом в качестве второго аргумента, можно с помощью функции
int MPI_Waitany(int count, MPI_Request array_of_ requests[], int *index, MPI_Status array_of_statuses[]);
индекс завершившейся операции заносится в аргумент index.
Ожидать завершения всех ранее запущенных операций пере-дач/приемов, квитанции на которые указаны в массиве, передаваемом в качестве второго аргумента, можно с помощью функции:
int MPI_Waitall(int count, MPI_Request array_of_ requestsf],
MPI_Status array_of_statuses[]).;
Ожидать завершения любого количества ранее запущенных операций передачи/приема, квитанции на которые указаны в массиве, передаваемом в качестве второго аргумента, можно с помощью функции
int MPI_Waitsome(int incount,
MPI_Request array_of_ requestsf],
94
int *outcount, int array_of_indices[], MPI_ Status array_of_statuses[]).;
Все функции MPI_Wa.it* и MPI_Test* автоматически освобождают квитанции завершившихся неблокирующих операций приема/передачи и устанавливают соответствующие дескрипторы в состояние MPI_REQUEST_NULL. Использовать такую квитанцию впоследствии нельзя, ее надо заново инициализировать.
Есть функция, с помощью которой можно узнать, завершилась ли операция приема/передачи без деструкции соответствующей квитанции:
int MPI_Request_get_status(MPI_Request request, int flag, MPI_Status *status);
Квитанция, связанная с любой неблокирующей операцией приема/передачи, может быть явно «освобождена» без ожидания завершения этой операции:
int MPI_Request_free(MPI_Request *request).;
Параметр request устанавливается в значение MPI_REQUEST_ NULL. Связанная с квитанцией операция не прерывается, однако проверить ее завершение с помощью MPI_Wa.it или MPI_ Test уже нельзя.
Любая блокирующая операция приема/передачи может быть инициирована в так называемом «отложенном» режиме. Это означает, что библиотека получает и запоминает аргументы этой операции, но начинает выполнять ее только при вызове одной из функций MPI_Start или MPI_Startall. Функции инициирования выглядят так:
int MPI_Send_init (void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm, MPI_Request *request).;
int MPI_Rsend_init (void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm, MPI_Request *request).;
int MPI_Ssend_init (void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm, MPI_Request *request).;
int MPI_Bsend_init (void *buf, int count, MPI_Datatype datatype,
int dest, int tag, MPI_Comm comm, MPI_Request *request);
95
int MPI_Recv_init (void *buf, int count, MPI_Datatype datatype,
int source, int tag, MPI_Comm comm, MPI_ Request *request);
Во всех функциях *send_init добавлен (по сравнению с соответствующей функцией *send) последний аргумент MPI_Request ★request. В функции MPI_Recv_init этот аргумент появился взамен аргумента MPI_Status ''"status.
Отложенные операции приема/передачи можно запустить на выполнение с помощью функций
int MPI_Start(MPI_Request *request);
int MPI_Startall(int count,
MPI_Request array_of_ requests[]);
В MPI есть группа функций, совмещающих операции приема и передачи. Они достаточно часто применяются при программировании «каскадных» или «линейных» схем, когда необходимо осуществлять однотипный обмен данными между ветвями. Примером является функция:
int MPI_Sendrecv (void* sendbuffer, int sendcount, MPI_Datatype senddatatype, int dest, int sendtag, void* recvbuffer, int recvcount, MPI_Datatype recvdatatype,
int src, int recvtag, MPI_Comm comm, MPI_Status* status);
Эта функция копирует данные из массива sendbuffer ветви с идентификатором src в буфер recvbuffer ветви с идентификатором dest.
Другая полезная функция:
int MPI_Sendrecv_replace (void* buffer, int count, MPI_Datatype datatype, int dest, int sendtag, int src, int recvtag, MPI_Comm comm, MPI_Status* status);
Данные также передаются из ветви src в ветвь dest, но используется только один буфер (вначале сообщение из него передается, потом в него же принимается).
6.4.	Коллективные операции взаимодействия процессов
Набор операций типа «точка-точка» в принципе является достаточным для программирования любых алгоритмов. Однако
96
во многих случаях бывает удобнее использовать так называемые «коллективные операции», которые, как правило, будут выполняться быстрее. Например, часто возникает потребность разослать значение переменной или массива из одного процессора всем остальным. Конечно, такую рассылку можно реализовать с использованием операций Send/Recv, однако гораздо удобнее воспользоваться специальной коллективной операцией.
Главное отличие коллективных операций от операций типа «точка-точка» состоит в том, что в них всегда участвуют все ветви указанного коммуникатора. Несоблюдение этого правила приводит либо к аварийному завершению программы, либо к еще более неприятному ее зависанию.
Отличительные особенности коллективных операций:
•	коллективные операции не взаимодействуют с операциями типа «точка-точка»;
•	коллективные операции всегда выполняются в режиме с блокировкой;
•	возврат из функции коллективного взаимодействия в каждой ветви происходит тогда, когда ее участие в коллективной операции завершилось, однако это не означает, что другие ветви завершили операцию;
•	количество получаемых данных должно быть в точности равно количеству посланных данных;
•	типы элементов посылаемых и получаемых сообщений должны совпадать;
•	сообщения не имеют идентификаторов (тегов).
Барьерная синхронизация всех ветвей коммуникатора:
int MPI_Barrier(MPI_Comm comm);
— исполнение программы продолжится только тогда, когда все ветви коммуникатора сотт вызовут эту функцию.
Рассылка одного и того же сообщения от одной ветви всем остальным ветвям данного коммуникатора:
int MPI_Bcast(void* buffer, int count, MPI_Datatype datatype, int sender, MPI_Comm comm);
Эта функция в ветви с номером sender работает как МР1_ Send, а во всех остальных ветвях — как MPI_Recv. Поэтому перед вызовом этой функции не надо выяснять с помощью условного оператора, является ветвь sender-ом или нет. После возврата
97
из этой функции массив buffer во всех ветвях будет содержать одинаковые значения в первых count элементах типа datatype.
Рассылка частей сообщений от одной ветви всем остальным ветвям данного коммуникатора:
int MPI_Scatter (void* sendbuffer, int sendcount,
MPI_Datatype senddatatype, void* recvbuffer, int recvcount, MPI_Datatype recvdatatype, int sender, MPI_Comm comm);
Эта функция тоже в ветви с номером sender работает как МР1_ Send, а во всех остальных ветвях—как MPI_Recv. Важно помнить, что должно выполняться условие recvcount^sizeof(recvdatatype') >= sendcount*sizeof(senddatatype'), иначе программа будет завершена аварийно. После возврата из этой функции массив recvbuffer будет содержать:
—	в ветви 0: первые sendcount (0 — sendcount-Г) элементов буфера sendbuffer передающей ветви;
—	в ветви 1: sendcount элементов с номерами от sendcount до 2*sendcount-l буфера sendbuffer передающей ветви.
Векторный вариант рассылки
int MPI_Scatterv(void *sendbuf, int *sendcnts, intdispls, MPI_Datatype senddatatype, void *recvbuf, int recvcnt, MPI_Datatype recvdatatype, int sender, MPI_Comm comm);
отличается от предыдущей функции тем, что количество элементов данных, рассылаемых по ветвям, определяется содержимым массива sendents (при этом для всех i должно быть recvcnt >= sendents [i], иначе параллельная программа аварийно завершается).
Сбор частей сообщений из всех ветвей коммуникатора в одну ветвь:
int MPI_Gather(void *sendbuf, int sendent,
MPI_Datatype senddatatype, void *recvbuf, int recvcnt, MPI_Datatype recvdatatype, int receiver, MPI_Comm comm);
Выполняемые действия в точности противоположны действиям функции MPI_Scatter. В буфер recvbuf ветви receiver складываются порции данных из всех остальных ветвей в строгом соответствии с их номерами. Значение аргумента recvcnt может быть больше или равно sendent, но не наоборот. При recvcnt < sendent программа будет завершена аварийно.
98
Векторный вариант сбора частей сообщений от всех ветвей данного коммуникатора в одну ветвь:
int MPI_Gatherv(void* sendbuf, int sendcount, MPI_Datatype senddatatype, void* recvbuf, int *recvcnts, int *displs, MPI_Datatype recvdatatype, int receiver, MPI_Comm comm);
Как обычно, для каждой принимаемой порции должно удовлетворяться условие: recvcnts[i]*sizeof(recvdatatype) >= sendcount'^sizeof(senddatatype), иначе — аварийное завершение.
Функция MPI_Allgather выполняется так же, как MPI Gather, но получателями являются все ветви группы данного коммуникатора. Данные, отправленные ветвью i из своего буфера sendbuf, помещаются в i-ю порцию буфера recvbuf каждой ветви. После завершения операции содержимое буферов приема recvbuf у всех ветвей будет одинаковым.
int MPI_Allgather(void* sendbuf, int sendcount, MPI_Datatype senddatatype, void* recvbuf, int recvcount, MPI_Datatype recvdatatype, MPI_Comm comm);
Функция MPI_Allgatherv является аналогом функции MPI_ Gatherv, но сборка данных выполняется всеми ветвями коммуникатора. Поэтому в списке аргументов отсутствует параметр receiver.
int MPI_Allgatherv(void* sendbuf, int sendcount, MPI_Datatype senddatatype, void* recvbuf, int *recvcnts, int *displs, MPI_Datatype recvdatatype, MPI_Comm comm);
Функция MPI_Alltoall совмещает в себе операции Scatter и Gather и является, по сути, расширением операции Allgather, когда каждая ветвь посылает различные данные разным получателям. Ветвь i посылает j-й блок своего буфера sendbuf ветви j, которая помещает его в i-й блок своего буфера recvbuf Количество посланных данных должно быть равно количеству полученных данных для каждой пары ветвей.
int MPI_Alltoall(void *sendbuf, int sendcount,
MPI_Datatype senddatatype, void *recvbuf, int recvcount, MPI_Datatype recvdatatype, MPI_Comm comm);
99
Функция MPI_Alltoallv является векторным вариантом функции MPI_Alltoall, позволяющим гибко управлять местоположением и размерами порций сообщений, передаваемых между ветвями.
int MPI_Alltoallv(void *sendbuf, int *sendcnts, int
*senddispls,
MPI_Datatype senddatatype, void *recvbuf, int *recvcnts, int *recvdispls, MPI_Datatype recvdatatype, MPI_Comm comm);
Еще более гибкой является функция MPI_Alltoallw, дополнительно позволяющая передавать/принимать данные разных типов:
int MPI_Alltoallw(void *sendbuf, int *sendcnts, int *senddispls, MPI_Datatype *senddatatypes, void *recvbuf, int *recvcnts, int *recvdispls, MPI_Datatype *recvdatatypes, MPI_Comm comm);
Операцией редукции называется операция, аргументом которой является вектор, а результатом — скалярная величина, полученная применением некоторой математической операции ко всем компонентам вектора. В частности, если компоненты вектора расположены в адресных пространствах ветвей, выполняющихся на различных процессорах, то в этом случае говорят о глобальной (параллельной) редукции.
Например, пусть в адресном пространстве ветвей некоторой группы имеются собственные копии переменной var (необязательно имеющие одно и то же значение). Тогда применение к этой переменной операции вычисления глобальной суммы или, другими словами, операции редукции SUM возвратит одно значение, которое будет содержать сумму всех локальных значений этой переменной.
Использование операций редукции является одним из основных средств организации распределенных вычислений. В функциях редукции MPI могут быть использованы предопределенные операции, перечисленные в табл. 6.2.
Таблица 6.2
Название	Операция	Разрешенные типы
МР1_МАХ	Максимум	С: int. FORTRAN: integer, floating point
MPI_MIN	Минимум	
100
Окончание табл. 6.2
Название	Операция	Разрешенные типы
MPI_SUM	Сумма	С: int. FORTRAN: integer, floating point, complex
MPI_PROD	Произведение	
MPI_LAND	Логическое AND	C int. FORTRAN: logical
MPI_LOR	Логическое OR	
MPI_LXOR	Логическое исключающее OR	
MPI_BAND	Поразрядное AND	C: int. FORTRAN: integer, byte
MPI_BOR	Поразрядное OR	
MPI_BXOR	Поразрядное исключающее OR	
MPI_ MAXLOC	Максимальное значение и его индекс	Специальные типы для этих функций
MPI_ MINLOC	Минимальное значение и его индекс	
Операции MAXLOC и MINLOC выполняются над специальными парными типами, каждый элемент которых хранит две величины: значения, по которым ищется максимум или минимум, и индекс элемента. В MPI имеется девять таких предопределенных типов, для языка С они перечислены в табл. 6.3.
Таблица 6.3.
Название	Типы
MPI_FLOAT_INT	float и int
MPI_DOUBLE_INT	double и int
MPI_LONG_INT	long и int
MPI_2INT	int и int
MPI_SHORT_INT	short и int
MPI_LONG_DOUBLE_INT	long double и int
Функция MPI_Reduce:
int MPI_Reduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, int receiver, MPI_Comm comm);
101
Над всеми элементами буфера sendbuf всех ветвей коммуникатора выполняется операция ор, результаты ее собираются в буфере recvbuf ветви receiver. В качестве операции ор можно использовать либо одну из предопределенных операций, либо операцию, сконструированную пользователем. Все предопределенные операции являются ассоциативными и коммутативными. Сконструированная пользователем операция, по крайней мере, должна быть ассоциативной. Порядок редукции определяется номерами ветвей в коммуникаторе. Тип datatype элементов должен быть совместим с операцией ор.
Похожа на MPI_Reduce функция MPI_Allreduce, только результаты собираются не в одной ветви, а во всех ветвях:
int MPI_AllReduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
Функция MPI-Reduces cutter отличается от MPI_Allreduce тем, что результат операции разрезается на непересекающиеся части по числу ветвей в группе, i-я часть посылается i-й ветви в ее буфер приема. Длины этих частей задает третий параметр, являющийся массивом.
int MPI_Reduce_scatter(void *sendbuf, void *recvbuf, int *recvcnts, MPI_Datatype datatype, MPI_Op op, MPI_ Comm comm)
Функция MPI_Scan выполняет префиксную включающую редукцию. Параметры такие же, как в MPI_Allreduce, но получаемые каждой ветвью результаты отличаются друг от друга. Операция пересылает в буфер приема i-й ветви редукцию значений из входных буферов ветвей с номерами 0, ..., i включительно.
int MPI_Scan(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
Функция MPI Exscan выполняет префиксную исключающую редукцию (собственные значения каждой ветви не участвуют в редукции):
int MPI_Exscan(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm);
102
Определение собственных операций для функций МР1_ Reduce, MPI_Reduce_scatter, MPI_Allreduce и MPIJScan:
int MPI_Op_create(MPI_User_function *function, int commute, MPI_Op *op);
Описание типа собственной функции выглядит следующим образом:
typedef void (MPI_User_function)(
void *a, void *b, int *len, MPI_Datatype *dtype);
Эта функция должна выполнять обработку данных из векторов а и b следующим образом:
b[i] =a[i] пользовательская операция> b[i]
для i=0, ..., len - 1.
Удаление ранее определенной собственной функции:
int MPI_Op_free(MPI_Op *ор);
В результате вызова переменная ор получает значение МР1_ OP_NULL.
6.5.	Производные типы и упаковка/распаковка данных
Все рассмотренные ранее коммуникационные операции позволяют посылать или получать последовательность элементов одного типа, занимающих смежные области памяти. При разработке параллельных программ часто возникает потребность передавать данные разных типов (например, структуры) или данные, расположенные в несмежных областях памяти ветви-отправителя (например, части массивов, не образующие непрерывную последовательность элементов).
Для эффективной пересылки данных в таких случаях MPI предоставляет два механизма:
•	возможность создания производных типов данных для использования в коммуникационных операциях вместо предопределенных типов MPI;
103
•	пересылку упакованных данных (процесс-отправитель упаковывает пересылаемые данные перед их отправкой, а процесс-получатель распаковывает их после получения).
Производные типы данных стандарта MPI не являются в полном смысле типами данных, как это понимается в языках программирования. Они не могут использоваться ни в каких других операциях, кроме операций передачи/приема сообщений. Производные типы данных MPI следует понимать просто как описатели расположения в памяти элементов базовых типов. Производный тип MPI представляет собой скрытый (opaque) объект, который специфицирует две вещи:
•	последовательность базовых типов;
•	последовательность их смещений.
Упорядоченный набор этих пар называется отображением (картой) типа:
Typemap = {(typeo, disp0),.., (type~, disp~ ~}.
Значения смещений в карте типа не обязательно должны быть неотрицательными, различными и упорядоченными по возрастанию.
Стандартный сценарий определения и использования производных типов включает следующие шаги:
•	производный тип строится из предопределенных типов MPI и ранее определенных производных типов с помощью функций-конструкторов;
•	новый производный тип регистрируется вызовом функции MPIJtype_commit;
•	после регистрации новый производный тип можно использовать в коммуникационных операциях и при конструировании других типов. Предопределенные типы MPI считаются зарегистрированными;
•	с производными типами можно выполнять некоторые дополнительные операции (в частности, узнавать их размер, протяженность и некоторые другие характеристики);
•	когда производный тип становится ненужным, он уничтожается функцией MPI_TypeJree.
Некоторые из перечисленных ниже функций-конструкторов появились только в стандарте MPI-2, однако для сохранения единства в изложении материала сведения о них приводятся здесь.
104
Самый простой конструктор типа MPI_Type_contiguous создает новый тип, элементы которого состоят из указанного числа элементов базового типа, занимающих смежные области памяти:
int MPI_Type_contiguous(int count, MPI_Datatype oldtype, MPI_Datatype *newtype);
Следующий конструктор создает тип, элемент которого представляет собой несколько равноудаленных друг от друга блоков из одинакового числа смежных элементов базового типа:
int MPI_Type_vector(int count, int blocklength, int stride, MPI_Datatype oldtype, MPI_Datatype *newtype);
Эта функция создает тип newtype, элемент которого состоит из count блоков, каждый из которых содержит одинаковое число blocklength элементов типа oldtype. Шаг stride между началом блока и началом следующего блока всюду одинаков и кратен протяженности представления базового типа.
Конструктор типа MPI_Type_hvector расширяет возможности конструктора MPI_Type_veclor, позволяя задавать произвольный шаг между началами блоков в байтах:
int MPI_Type_create_hvector(int count, int blocklength, MPI_Aint stride, MPI_Datatype oldtype, MPI_Datatype *newtype);
Конструктор типа MPI_Type-indexed является более универсальным по сравнению с MPI_Type_veclor, так как элементы создаваемого типа состоят из произвольных по длине блоков с произвольным смещением блоков от начала размещения элемента. Смещения задаются в элементах базового типа:
int MPI_Type_create_indexed(int count,
int *array_of_ blocklengths,
int *array_of_displacements,
MPI_Datatype oldtype, MPI_Datatype *newtype);
Эта функция создает тип newtype, каждый элемент которого состоит из count блоков, где i-й блок содержит array_of_ blocklengths[i] элементов базового типа и смещен от начала размещения элемента нового типа на array_of-displacements [i] элементов базового типа.
105
Конструктор типа MPI_Type_create_hindexed идентичен конструктору MPI_Type-indexed за исключением того, что смещения измеряются в байтах:
int MPI_create_hindexed(int count,
int *array_of_ blocklengths,
MPI_Aint *array_of_displacements, MPI_ Datatype oldtype, MPI_Datatype *newtype).;
Элемент нового типа состоит из count блоков, где i-й блок содержит array_of-blocklengths [Г\ элементов старого типа и смещен от начала размещения элемента нового типа на array_of_ displacements[i] байт.
Конструктор типа MPI_Type_create_indexed_block похож на конструктор MPI_Type-indexed, за исключением того, что все блоки одинаковы:
int MPI_Type_create_indexed_block(int count,
int blocklength,
MPI_Aint *array_of_displacements, MPI_ Datatype oldtype, MPI_Datatype *newtype);
Конструктор типа MPI_Type_create_struct — самый универсальный из всех конструкторов. Создаваемый им тип является структурой, состоящей из произвольного числа блоков, каждый из которых может содержать произвольное число элементов одного из базовых типов и может быть смещен на произвольное число байтов от начала размещения структуры:
int MPI_Type_create_struct(int count,
int *array_of_ blocklengths,
MPI_Aint *array_of_displacements,
MPI_ Datatype *array_of_types, MPI_Datatype *newtype);
Эта функция создает тип newtype, элемент которого состоит из count блоков, где i-й блок содержит array_of-blocklengths [i] элементов типа array_of_types[i]. Смещение i-ro блока от начала размещения элемента нового типа измеряется в байтах и задается в array_of_displacements[i].
Конструктор типа данных «субмассив»:
int MPI_Type_create_subarray(int ndims,
int* array_ of_sizes,
106
int* array_of_subsizes, int* array_of_starts, int order, MPI_Datatype oldtype, MPI_Datatype *newtype);
создает тип данных MPI, описывающий n-мерный субмассив n-мерного массива. Субмассив может находиться в любом месте полного массива и может быть любого ненулевого размера вплоть до размера полного массива. Этот конструктор позволяет создавать типы файлов для доступа к массивам, разбитым между процессами по блокам через один файл, содержащий глобальный массив.
И, наконец, конструктор распределенного массива поддерживает распределения данных, сходные с HPF (High Performance Fortran). Кроме этого, в отличие от HPF, порядок хранения может быть задан как для массивов Си, так и для ФОРТРАНа:
int MPI_Type_create_darray(int size,
int rank, int ndims,
int array_of_gsizes[], int array_of_distribs[],
int array_of_dargs[], int array_of_psizes[], int order, MPI_Datatype oldtype, MPI_Datatype* newtype):
Функция MPI_Type_commit регистрирует созданный производный тип. Только после регистрации новый тип может использоваться в коммуникационных операциях:
int MPI_Type_commit(MPI_Datatype *datatype).;
Любой зарегистрированный тип данных можно продублировать с помощью функции
int MPI_Type_dup(MPI_Datatype datatype, MPI_Datatype *newtype);
Функция
int MPI_Type_free(MPI_Datatype *datatype).;
уничтожает производный тип. Эта функция устанавливает дескриптор типа в состояние MPI_DATATYPE_NULL, что не повлияет ни на выполняющиеся в данный момент коммуникационные операции с этим типом данных, ни на производные типы, которые ранее были определены через уничтоженный тип.
107
Передача/прием упакованных данных
Функция МР1_Раск упаковывает элементы предопределенного или производного типа, помещая их побайтное представление в выходной буфер:
int MPI_Pack(void* inbuf, int incount,
MPI_Datatype datatype,
void *outbuf, int outsize, int *position, MPI_Comm comm);
Функция MPI_Pack упаковывает incount элементов типа datatype из области памяти с начальным адресом inbuf. Результат упаковки помещается в выходной буфер с начальным адресом outbufn размером outsize байт. Параметр position указывает текущую позицию в байтах, начиная с которой будут размещаться упакованные данные. После возврата из функции значение position будет увеличено на число упакованных байт, указывая на первый свободный байт в выходном буфере. Параметр сотт при последующей посылке упакованного сообщения будет использован как коммуникатор.
Функция MPI_Pack_size позволяет определить размер буфера, необходимый для упаковки заданного количества данных типа datatype:
int MPI_Pack_size(int incount, MPI_Datatype datatype, MPI_Comm comm, int *size);
Функция MPI_Unpack извлекает заданное число элементов некоторого типа из побайтного представления элементов во входном массиве:
int MPI_Unpack(void* inbuf, int insize, int *position, void *outbuf, int outcount, MPI_Datatype datatype, MPI_Comm comm);
Эта функция извлекает outcount элементов типа datatype из побайтного представления элементов в массиве inbuf, начиная с адреса position, и помещает их в область памяти с начальным адресом outbuf. После возврата из функции параметр position будет увеличен на размер распакованного сообщения.
б.б.	Управление группами ветвей и коммуникаторами
Коммуникаторы определяют области действия любых операций обмена данными в MPI.
108
Типичной проблемой, которую может решить использование коммуникаторов, является проблема недопущения «пересечения» обменов по тегам. В самом деле, если программист использует какую-то библиотеку параллельных методов, то он часто не знает, какие теги использует эта библиотека при передаче сообщений. В таком случае существует опасность, что тег, выбранный программистом, совпадет с одним из тегов, используемых библиотекой.
Чтобы этого не произошло, программист может создать собственный коммуникатор и выполнять операции обмена в его рамках. Несмотря на то что этот коммуникатор будет содержать те же ветви программы и такие же теги, что и коммуникатор библиотеки, сообщения, передаваемые в рамках этих двух разных коммуникаторов, не перепутаются.
Управление группами
С понятием коммуникатора тесно связано понятие группы ветвей. Под группой понимают упорядоченное множество ветвей программы. Каждой ветви в группе соответствует уникальный номер — ранг. Группа — отдельное понятие MPI, и операции с группами могут выполняться отдельно от операций с коммуникаторами, но операции обмена для указания области действия всегда используют коммуникаторы, а не группы. С группами ветвей в MPI допустимы следующие действия:
•	объединение групп;
•	пересечение групп;
•	образование разности групп.
Новая группа может быть создана только из уже существующих групп. В качестве исходной группы при создании новой может быть использована, например, группа, связанная с предопределенным коммуникатором MPI_COMM_WORLD. При ручном конструировании групп может оказаться полезной специальная пустая группа MPI_COMM_EMPTY.
Для доступа к группе ветвей, связанной с коммуникатором comm, используется функция.
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group);
Размер группы (количество ветвей в ней) можно получить с помощью функции.
int MPI_Group_size(MPI_Group group, int *size);
109
Номер ветви в группе можно определить, вызвав функцию int MPI_Group_rank(MPI_Group group, int *rank);
Если есть две группы ветвей (неважно, каким образом они были образованы), то с помощью функции
int MPI_Group_translate_ranks(MPI_Group groupl,
int n, int *ranksl,
MPI_Group group2, int *ranks2);
можно установить соответствие между номерами ветвей в этих группах.
Две группы ветвей можно сравнить:
int MPI_Group_compare(MPI_Group groupl, MPI_Group group2, int *result).;
где result — переменная, получающая в результате сравнения значение:
—	MPI_IDENT, если все элементы двух групп одинаковы и следуют друг за другом в одинаковом порядке;
—	MPIJSIMILAR, если все элементы двух групп одинаковы, но порядок следования различен;
—	MPIJUNEQUAL — во всех остальных случаях.
Однажды созданная группа не может быть изменена. Она может быть только уничтожена (освобождена). Для конструирования новых групп можно использовать следующие функции.
Создание новой группы из явно указанных ветвей старой группы:
int MPI_Group_incl(MPI_Group oldgroup, int n, int *ranks, MPI_Group *newgroup);
Эта функция создает новую группу newgroup, которая состоит из ветвей существующей группы, перечисленных в массиве ranks. Ветвь с номером i в новой группе есть ветвь с номером ranks[z] в существующей группе. Каждый элемент в массиве ranks должен иметь корректный номер ветви из группы oldgroup, и среди этих элементов не должно быть совпадающих, иначе вся программа завершится аварийно.
110
Создание новой группы из таких ветвей старой группы, которые не принадлежат явно указанному перечню:
int MPI_Group_excl(MPI_Group oldgroup, int n, int *ranks, MPI_Group *newgroup);
Эта функция создает новую группу newgroup, которая состоит из всех ветвей существующей группы, за исключением ветвей, перечисленных в массиве ranks. Каждый элемент в массиве ranks должен иметь корректный номер ветви из группы oldgroup, и среди этих элементов не должно быть совпадающих, иначе вся программа завершится аварийно.
Есть две функции, являющиеся обобщением функций МР1_ Group_incl и MPI_Group_excl:
int MPI_Group_range_incl(MPI_Group oldgroup, int n,
int ranges[][3], MPI_Group *newgroup);
int MPI_Group_range_excl(MPI_Group oldgroup, int n,
int ranges[][3],
MPI_Group *newgroup);
У этих функций одномерные массивы номеров ranks заменены двумерными массивами триплетов, каждый из которых имеет вид:
<нижняя граница> <верхняя граница> <шаг>
Границы и шаг должны определять допустимые номера ветвей, причем все номера ветвей новой группы, определяемые с их использованием, должны быть различными, иначе вся программа завершится аварийно. Последовательность, определяемая триплетом, может быть пустой (например, если верхняя граница меньше нижней). Эти функции удобно применять в программах для суперкомпьютеров с очень большим (десятки и сотни тысяч) количеством ядер. Пример массива ranges:
int rngs[][] = {{2, 32768, 16},{0, 131072, 256}};
Здесь имеются два триплета:
1) 2, 32768,16 — порождает последовательность номеров 2, 18, 34, 50, 66, ..., 32 738, 32 754 (всего 2047 номеров);
2) 0, 131072, 256 — порождает последовательность 0, 256, 512, ..., 130 816, 131 072 (всего 512 номеров).
ill
Следующие три функции имеют одинаковые наборы аргументов и создают новую группу как результат заданной теоретико-множественной операции над множествами ветвей двух групп.
int MPI_Group_union(MPI_Group groupl, MPI_Group group2, MPI_Group *newgroup);
int MPI_Group_intersection(MPI_Group groupl,
MPI_ Group group2, MPI_Group *newgroup);
int MPI_Group_difference(MPI_Group groupl, MPI_Group group2, MPI_Group *newgroup);
Уничтожение группы выполняется путем вызова функции
int MPI_Group_free(MPI_Group *group)
Управление коммуникаторами
Создание коммуникатора — операция коллективная (и должна вызываться всеми процессами коммуникатора).
Функция MPI_Comm_dup копирует уже существующий коммуникатор:
int MPI_Comm_dup(MPI_Comm oldcomm, MPI_comm *newcomm);
В результате вызова этой функции будет создан новый коммуникатор newcomm, включающий в себя те же ветви, что и коммуникатор oldcomm.
Для создания нового коммуникатора служит функция МР1_ Сотт create:
int MPI_comm_create(MPI_Comm oldcomm, MPI_Group group, MPI_Comm *newcomm);
Этот вызов создает новый коммуникатор newcomm, который будет включать в себя ветви группы group коммуникатора oldcomm. Это коллективная операция, и она должна вызываться всеми ветвями родительского коммуникатора, даже если ветвь не входит в группу нового коммуникатора. Такие ветви в качестве нового коммуникатора получат MPI_COMM_NULL.
Функция расщепления коммуникатора:
int MPI_Comm_split(MPI_Comm oldcomm, int color, int key, MPI_Comm *newcomm).;
112
Эта функция расщепляет группу, связанную с родительским коммуникатором, на непересекающиеся подгруппы, по одной на каждое значение признака подгруппы color. Значение color должно быть неотрицательным. Каждая подгруппа будет содержать ветви с одним и тем же значением признака color. Параметр key управляет упорядочиванием внутри новых групп: меньшему значению key соответствует меньшее значение идентификатора ветви. В случае равенства параметра key для нескольких ветвей упорядочивание выполняется в соответствии с порядком в родительской группе. Каждая ветвь получит свой новый коммуникатор.
Функция сравнения двух коммуникаторов:
int MPI_Comm_compare(MPI_Comm comml, MPI_Comm comm2, int *result).;
Функция уничтожения коммуникатора:
int MPI_Comm_free(MPI_Comm *comm);
6.7.	Управление виртуальными топологиями
Топология — механизм MPI, который позволяет устанавливать дополнительную систему адресации для процессов. Топология в MPI является логической и никак не связана с топологией физической среды передачи данных. Введение в библиотеку MPI поддержки топологий связано с тем, что большинство прикладных алгоритмов устроено таким образом, что ветви должны быть упорядоченными в соответствии с некоторой топологией. В MPI поддерживаются два вида топологий — прямоугольная решетка произвольной размерности (декартова топология) и граф.
Декартовы топологии или решетки часто применяются при решении прикладных задач. Известно много алгоритмов, в которых используются «сетки». Один из наиболее широко представленных классов таких задач — сеточные методы решения дифференциальных уравнений в частных производных. Именно для упрощения программирования «сеточных» задач в MPI и были введены функции управления топологиями.
113
Чтобы создать топологию вида «решетка», нужно воспользоваться функцией:
int MPI_Cart_create(MPI_Comm oldcomm, int ndims, int *dims, int *periods, int reorder, MPI_Comm grid_comm);
Операция создания топологии является коллективной операцией, ее должны выполнить все процессы коммуникатора oldcomm. Если какие-то ветви при этом не попадают в новую группу, для них возвращается результат MPI_COMM_NULL. В случае, когда размеры заказываемой сетки больше имеющегося в группе числа ветвей, функция завершается аварийно. Значение параметра reorder = false означает, что идентификаторы всех ветвей в новой группе будут такими же, как в старой группе, при значении true MPI будет пытаться перенумеровать их с целью оптимизации коммуникаций. С помощью этой функции можно создавать топологии с произвольным числом измерений, причем по каждому измерению в отдельности можно накладывать различные граничные условия. Таким образом, для одномерной топологии мы можем получить или линейную структуру, или кольцо, в зависимости от того, какие граничные условия будут наложены. Для двумерной топологии — соответственно либо прямоугольник, либо цилиндр, либо тор и т. д.
Перед вызовом функции MPI_Cart_create полезно определить оптимальное распределение ветвей программы по предполагаемой решетке с помощью вызова функции
int MPI_Dims_create(int nnodes, int ndims, int *dims).;
перед вызовом функции в массив dims должны быть занесены целые неотрицательные числа. Если элементу массива dims [i] присвоено положительное число, то для этой размерности решетки вычисление не производится (число ветвей вдоль этого направления считается заданным). Вычисляются только те компоненты dims [i], для которых перед обращением к функции были присвоены значения 0. Функция стремится создать максимально равномерное распределение ветвей вдоль направлений, выстраивая их по убыванию.
Для определения того, связана ли (и какая) топология с коммуникатором, служит функция
int MPI_Topo_test(MPI_Comm comm, int *status);
114
Функция MPI_Topo_test возвращает через переменную status топологию коммуникатора сотт. Возможные значения:
•	MPI_CART — декартова топология;
•	MPI_GRAPH — топология графа;
•	MPI_UNDEFINED — топология не задана.
Функция опроса числа измерений декартовой топологии:
int MPI_Cartdim_get(MPI_Comm comm, int *ndims).;
Результат, возвращаемый этой функцией, может быть использован в качестве параметра для вызова функции MPI_Cart_get, которая служит для получения более детальной информации.
int MPI_Cart_get(MPI_Comm comm, int ndims, int *dims, int *periods, int *coords).;
Для определения декартовых координат ветви по ее рангу (номеру) можно воспользоваться функцией
int MPI_Cart_coords(MPI_Comm comm,int rank,int maxdims, int *coords).;
Функция MPI_Cart_rank возвращает ранг ветви по ее декартовым координатам:
int MPI_Cart_rank(MPI_Comm comm, int *coords, int *rank);
Для совокупностей ветвей, организованных в гиперкуб, могут выполняться обмены особого рода — сдвиги. Существует два вида сдвигов:
1)	циклический сдвиг на к элементов вдоль ребра решетки. Данные от ветви с номером п пересылаются ветви с номером (п + к) % dim, где dim — размерность измерения, вдоль которого производится сдвиг;
2)	линейный сдвиг на к позиций вдоль ребра решетки. Данные от ветви с номером п пересылаются ветви с номером (и + к) (если такая существует).
Собственно передачи/приемы данных выполняются функцией MPI_Sendrecv, но вычисление номеров ветвей приемников/по-лучателей целесообразно выполнять с помощью вызова функции
int MPI_Cart_shift(MPI_Comm comm, int direction, int displ, int* source, int* dst);
115
Часто используемая операция — выделение в декартовой топологии подпространств меньшей размерности и связывание с ними отдельных коммуникаторов. Функция выделения подпространства в декартовой топологии:
int MPI_Cart_sub(MPI_Comm comm, int *remain_dims, MPI_Comm *newcomm);
Создание коммуникаторов с произвольной топологией выполняется функцией
int MPI_Graph_rreate(MPI_Comm oldcomm, int nnodes,
int *index, int *edges, int reorder, MPI_Comm *comm_graph);
Эта функция передает дескриптор новому коммуникатору, к которому присоединяется информация о графовой топологии. Если размер декартовой решетки nnodes меньше, чем размер группы коммуникатора, то некоторым ветвям возвращается значение MPI_COMM_NULL, по аналогии с MPI_Cart_spl.it и MPI_Comm_split. Вызов будет неверным, если он определяет граф большего размера, чем размер группы исходного коммуникатора.
Структуру графа определяют три параметра: nnodes, index и edges. Способ определения аргументов nnodes, index и edges иллюстрируется следующим простым примером.
Пусть есть четыре ветви, имеющие следующих соседей: ветвь 0 — 1 и 3, ветвь 1 — 0, ветвь 2 — 3 и ветвь 3 — 0 и 2. Тогда исходными аргументами функции должны быть:
nnodes = 4, index = {2, 3, 4, 6}, edges = {1, 3, 0, 3, 0, 2}
Заметим, что не имеет никакого смысла создавать виртуальную топологию полного графа. Возможность передавать данные между любой парой ветвей программы возникает изначально, как только разработчик подключает к программе библиотеку MPI. Виртуальные топологии — это средство упростить формирование координат «соседей» ветви для указания их в вызовах функций приема/передачи данных и не более того.
Для коммуникаторов с графовой топологией есть функции получения сведений, аналогичные функциям MPI_Cartdims_get и MPI_Cart_get. Размеры массивов вершин и ребер графа могут быть получены с помощью вызова функции
116
int MPI_Graphdims_get(MPI_Comm comm, int *nnodes, int *nedges);
Информация, получаемая в результате вызова этой функции, может быть использована для корректного определения размера векторов index и edges в последующем вызове функции
int MPI_Graph_get(MPI_Comm comm, int nnodes, int nedges, int *index, int *edges).;
Информацию о смежных ветвях любой вершины графа можно получить с помощью двух функций. Первая из них является вспомогательной и возвращает количество «соседей» ветви с заданным рангом:
int MPI_Graph_neighbors_count(MPI_Comm comm, int rank, int *nneighbors).;
После получения количества соседей можно получить полный перечень этих соседей с помощью функции:
int MPI_Graph_neighbors(MPI_Comm comm, int rank, int maxneighbors, int *neighbors);
При необходимости использования других функций можно использовать руководства по программированию для MPI и справочные материалы.
6.8.	Удаленный доступ к памяти (односторонние взаимодействия)
Удаленный доступ к памяти (RMA) расширяет механизмы взаимодействий MPI, позволяя одной ветви определять все коммуникационные параметры как для посылающей стороны, так и для получающей.
Во время любого из рассмотренных ранее коллективных MPI-взаимодействий взаимосвязанно выполняются:
•	передача данных от отправителя к получателю;
•	синхронизация отправителя и получателя.
В отличие от коллективных взаимодействий, функции RMA спроектированы таким образом, что эти операции разделяются, передача данных в память / из памяти одной ветви выпол
117
няется без ее участия по запросу другой ветви. Именно поэтому такие передачи и называются односторонними взаимодействиями.
Инициатором RMA-взаимодействия называется ветвь, которая запрашивает доступ к памяти другой ветви, и адресатом — ветвь, к памяти которой выполняется обращение. Инициатор и адресат могут быть одной и той же ветвью, тогда операции передачи данных используются просто для их перемещения в пределах памяти этой ветви.
Тем не менее коллективность имеется и в односторонних взаимодействиях. Прежде всего коллективными являются: операция создания «окна» для удаленного доступа к памяти некоторой ветви параллельной программы, операция удаления (уничтожения) такого окна и операции синхронизации.
Собственно удаленный доступ к окну памяти реализуется с использованием одной из трех функций:
1)	MPI_Put передает данные из памяти инициатора в память адресата;
2)	MPI_Get передает данные из памяти адресата в память инициатора;
3)	MPI_Accumulate модифицирует данные в памяти адресата (например, добавляя к ним значения, посланные из памяти инициатора).
Эти три операции являются неблокирующими: т. е. вызов функции инициирует взаимодействие, после чего передача данных может продолжаться после возврата из вызова. Выполнение передачи завершается как в инициаторе, так и в адресате, когда инициатором выдан последующий синхронизационный вызов к участвующему оконному объекту.
Коллективная операция инициализации позволяет каждой ветви из группы интракоммуникатора определить «окно» в своей памяти, которое становится доступным для других ветвей этой группы. Вызов возвращает объект со скрытой структурой, представляющий группу процессов, которые являются обладателями и имеют доступ к окну и к его атрибутам.
int MPI_Win_create(void *base, MPI_Aint size, int disp_unit, MPI_Info info, MPI_Comm comm, MPI_Win *win);
Этот вызов возвращает оконный объект, который может использоваться всеми ветвями коммуникатора сотт для выпол
118
нения RMA-операций. Каждая ветвь определяет окно в своей памяти, которое она предоставляет для дистанционного доступа всем ветвям из группы сотт. Окно состоит из size байт, начинающихся с адреса base. Если значение size равно нулю, то окно из памяти этой ветви недоступно другим ветвям. Для упрощения адресной арифметики в будущих RMA-операциях задается аргумент disp_unit, определяющий коэффициент масштабирования. Обычным выбором для значения disp_unit являются либо 1 (нет масштабирования), либо sizeof(type) для окон, которые состоят из массива элементов типа type. В последнем случае в RMA-вызовах можно будет использовать индексы элементов массивов, причем правильное масштабирование к байтовому смещению обеспечивается даже в неоднородной среде. Параметр info предоставляет подсказку исполняющей системе для выполнения оптимизации относительно ожидаемой структуры использования окна. Разные ветви в группе коммуникационного взаимодействия сотт могут определять окна-адресаты, полностью различающиеся по расположению, размеру, по единицам смещения и параметру (аргументу) info. Одна и та же область памяти может использоваться одновременно в нескольких разных окнах. Однако следует помнить, что перекрывающиеся окна являются источником трудно обнаруживаемых ошибок. Окно можно создать в любой части памяти ветви. На некоторых платформах эффективность окон будет выше, если они создаются в памяти, которая распределена вызовом МР1_А11ос_тет. Эффективность также обычно улучшается, когда границы окна выравниваются по «естественным» границам (слово, двойное слово, строка кэша, страничный блок и т. д.).
Ранее созданное окно может быть уничтожено с помощью коллективной операции:
int MPI_Win_free(MPI_Win *win);
Эта операция освобождает оконный объект и заносит в переменную win пустой дескриптор (значение MPI_WIN_NULL). Эта функция может вызываться ветвью только после того, как она завершила свое участие в RMA-взаимодействиях с оконным объектом win:
—	ветвь вызвала MPI-WinJvnce или MPI_Win_wait, чтобы выполнить согласование с предыдущим вызовом MPI_Win_post;
—	ветвь вызвала MPI_Win_complete, чтобы выполнить согласование с предыдущим вызовом MPI_Win_start;
119
— ветвь выполнила вызов MPI_Win_unlock, чтобы выполнить согласование с предыдущим вызовом MPI_Win_lock.
Функция MPI_WinJree внутри использует барьерную синхронизацию: никакая из ветвей не получит возврат из нее, пока все ветви из группы данного окна win не вызовут эту функцию.
Операция передачи данных из памяти инициатора в память адресата:
int MPI_Put(void *origin_addr, int origin_count, MPI_Datatype origin_datatype, int target_rank, MPI_Aint target_disp, int target_count, MPI_Datatype target_ datatype, MPI_Win win);
Из памяти ветви-инициатора, начиная с адреса origin_addr, передаются origin_count следующих друг за другом записей типа origin_datatype, в память адресата, определяемого парой win и target_rank. Данные записываются в буфер адресата по адресу target_addr, равному winjbase + target_disp * disp_ unit, где win_base и disp_unit — базовый адрес и единица смещения, определенные при инициализации оконного объекта win ветвью-адресатом. Передача данных происходит так же, как если бы инициатор выполнил операцию MPIJSend с аргументами origin_addr, origin_count, origin_datatype, target rank, tag, comm, и процесс-получатель выполнил операцию MPI_Recv с аргументами target_addr, tar get-datatype, source, tag, comm, где targetjaddr — это адрес буфера получателя, вычисленный, как указано в предыдущем абзаце, а сотт — это коммуникатор для группы ветвей win.
Операция передачи данных из памяти адресата в память инициатора:
int MPI_Get(void *origin_addr, int origin_count, MPI_Datatype origin_datatype, int target_rank, MPI_Aint target_disp, int target_count, MPI_Datatype target_ datatype, MPI_Win win);
Все аргументы точно такие же, как в предыдущем случае. Передача данных выполняется в обратном направлении: из памяти адресата в память инициатора.
Модификация данных в памяти адресата с использованием заданной операции и данных из памяти инициатора:
120
int MPI_Accumulate(void *origin_addr, int origin_ count, MPI_Datatype origin_datatype, int target_rank, MPI_Aint target_disp, int target_count, MPI_Datatype target_datatype, MPI_Op op, MPI_Win win);
Все аргументы такие же, как и у двух предыдущих функций, добавлен только аргумент ор, используемый точно так же, как в операции MPI_Reduce для «свертки» значений данных из буферов ветви-инициатора и ветви-адресата.
Синхронизация при удаленном доступе к памяти
RMA-взаимодействия могут выполняться в двух режимах:
1) с активным адресатом, когда данные перемещаются из памяти одной ветви в память другой, и обе явно участвуют в коммуникациях. Такая модель похожа на обычную передачу сообщений, за исключением того, что все аргументы передачи данных обеспечиваются одной ветвью, а вторая участвует только в синхронизации;
2) с пассивным адресатом, когда данные перемещаются из памяти одной ветви в память другой, и только инициатор явно участвует в передаче. При этом два или более инициатора могут взаимодействовать, имея доступ к одним и тем же данным в окне ветви-адресата. Ветвь, владеющая адресным окном, может отличаться от взаимодействующих ветвей, в этом случае она явно не участвует в коммуникациях. Этот режим близок к модели общей (shared) памяти, где общие данные могут быть доступны всем ветвям, невзирая на их расположение.
И в том, и в другом режиме должна обеспечиваться синхронизация действий ветвей, обеспечивающая единственность периода доступа к любой порции данных в любом окне. Под периодом доступа какой-либо ветви к какому-либо буферу памяти понимается промежуток времени, в течение которого только эта ветвь имеет право читать/изменять данные в этом буфере. Нарушение этого принципа обычно влечет за собой трудно обнаруживаемые ошибки в программе. MPI обеспечивает три механизма синхронизации.
1.	Коллективная синхронизация MPI-WinJence. Она обеспечивает простую модель синхронизации, которая часто используется при параллельных вычислениях: а именно слабосинхронную модель (loosely synchronous model), при которой вычислительные фазы перемежаются с фазами коммуникаций. Такой механизм более всего пригоден для слабосинхронных
121
алгоритмов, когда граф взаимодействующих процессов меняется очень часто либо когда каждый процесс взаимодействует со многими другими. Этот вид синхронизации используется для коммуникаций с активным адресатом. Периоды доступа начинаются и заканчиваются вызовами функции MPI-WinJence.
2.	Четыре функции, MPI_Win_start, MPI_Win_complete, MPI_ Win_post и MPI_Win_wait, могут использоваться, чтобы свести затраты на синхронизацию к минимуму. Синхронизируются только пары взаимодействующих процессов, и это происходит только тогда, когда синхронизация необходима, чтобы корректно упорядочить RMA-обращения к окну (с учетом локальных обращений к этому же окну).
3.	И последний механизм, блокировки, обеспечиваются двумя функциями: MPI_Win_lock и MPI_Win_unlock. Синхронизация с блокировками полезна для таких приложений, которые эмулируют модель с общей памятью через MPI-вызовы. Типичный пример — модель «доска объявлений», где ветви программы могут читать или обновлять различные части «доски объявлений» в случайные промежутки времени.
Функция коллективной «заборной» синхронизации:
int MPI_Win_fence(int assert, MPI_Win win);
Здесь assert задает контекст вызова для возможной внутренней оптимизации коммуникационных процессов MPI путем комбинирования флагов (допустимо отсутствие флагов, т. е. значение 0):
•	MPI_MODE_NOSTORE — локальная ветвь не меняла данные окна с момента последней синхронизации. Отменяет необходимость синхронизации кэшей;
•	MPI_MODE_NOPUT — данные локального окна не будут изменены вызовами MPI_Put и MPI_Accumulate до вызова МР1_ Win_wait; позволяет избежать синхронизации кэшей в момент вызова MPI_Win_wait;
•	MPI_MODE_NOPRECEDE — до этого вызова не было других вызовов RMA (имеется в виду — после последнего вызова MPI-WinJbnce). Если одна из ветвей окна устанавливает этот флаг, все остальные ветви группы должны тоже установить его;
•	MPI_MODE_NOSUCCEED — после этого вызова не будет других вызовов RMA (до следующей коллективной синхронизации). Если одна из ветвей окна устанавливает этот флаг, то все остальные должны тоже установить его.
122
Вызов функции MPI-WmJbnce подразумевает выполнение барьерной синхронизации всех ветвей, входящих в группу окна win.
Общая синхронизация с активным адресатом. Начало периода доступа к окну из группы ветвей-инициаторов:
int MPI_Win_start(MPI_Group group, int assert, MPI_ Win win);
Здесь assert задает контекст вызова для возможной внутренней оптимизации коммуникационных процессов MPI:
•	MPI_MODE_NOCHECK — соответствующий вызов МР1_ Win_post уже завершился на всех удаленных процессах до момента вызова MPI_Win_start. Этот флаг может быть указан в MPI_Win_start тогда и только тогда, когда он указан в соответствующем вызове MPI_Win_post;
•	значение assert = 0 может быть использовано в любой паре связанных вызовов MPI_Win_start и MPI_Win_post.
Вызов этой функции начинает период RMA-доступа ветвей группы group к данным окна win с активным адресатом. Каждый адресат должен вызвать функцию MPI_Win_post с соответствующими (такими же) аргументами. До тех пор, пока это не будет сделано, любое RMA-обращение к окну-адресату будет задержано.
Завершение периода доступа к окну активного адресата:
int MPI_Win_complete(MPI_Win win);
Вызов этой функции завершает период RMA-доступа к окну win, начатый вызовом MPI_Win_start. Все коммуникационные RMA-вызовы к окну win, созданные во время этого периода, завершатся в инициаторе к моменту, когда произойдет возврат из вызова этой функции. MPI_Win_complete заставляет завершиться предшествующие RMA-вызовы в инициаторе, но не в адресате. Вызов MPI_put или MPI_accumulate может еще не выполниться у адресата, в то время как он уже выполнился у инициатора.
Общая синхронизация с активным адресатом. Начало периода доступа к локальному окну в ветви-адресате:
int MPI_Win_post(MPI_Group group, int assert, MPI_Win win);
Здесь assert задает контекст вызова для возможной внутренней оптимизации коммуникационных процессов MPI:
123
•	MPI_MODE_NOCHECK — соответствующие функции MPI_ Win_ start еще не были вызваны ни одной из ветвей, намеревающейся осуществить доступ к данным до момента вызова MPI_Win_post;
•	MPI_MODE_NOSTORE — локальная ветвь не изменила данные окна с момента последней синхронизации. Отменяет необходимость синхронизации кэшей;
•	MPI_MODE_NOPUT — данные локального окна не будут изменены вызовами MPI_Put и MPI_Accumulate до вызова МР1_ Win_wait.
Начинает период предоставления RMA-доступа для локального окна, связанного с win. Только ветви из группы group будут иметь доступ к данным окна при помощи RMA-вызовов во время этого периода. Каждая ветвь в группе должна обеспечить соответствующий вызов MPI_Win_start. Функция MPI_Win_post не блокирует вызывающую ее ветвь программы.
Завершение периода доступа к локальному окну:
int MPI_Win_wait(MPI_Win win);
завершает период предоставления RMA-доступа к окну win, начатый вызовом MPI_Win_post. Этот вызов соответствует вызовам MPI_Win_complete, созданным каждым инициатором, которые имели доступ к окну. Вызов MPI_Win_wait будет блокировать ветвь, пока не завершатся все соответствующие вызовы MPI_Win_complete. Это гарантирует, что все инициаторы закончили свой RMA-доступ к локальному окну.
Неблокирующая версия функции MPI_Win_wait:
int MPI_Win_test(MPI_Win win, int flag);
Этот вызов является неблокирующей версией MPI_Win_wait. Он возвращает flag = true, если из вызова MPI_Win_wail может быть выполнен возврат, и flag = false в противном случае. Эффект от возвращения MPI_Win_test с flag = true такой же, как эффект от возвращения MPI_Win_wait. Если же возвращен flag = =false, тогда у вызова нет видимого эффекта. MPI_Win_test должен вызываться только там, где можно вызвать MPI_Win_wait. Как только произойдет возврат с кодом flag = true для некоторого оконного объекта, MPI_Win_test больше не должен вызываться для этого окна, пока оно не будет снова предоставлено для доступа.
124
Захват окна одной ветвью-инициатором:
int MPI_Win_lock(int lock_type, int rank, int assert, MPI_Win win);
Здесь lockjtype — тип захвата:
•	MPI_LOCK_EXCLUSIVE — только для ветви с номером rank;
•	MPI_LOCK_SHARED — возможен доступ из любой ветви;
•	assert — задает контекст вызова для возможной внутренней оптимизации коммуникационных процессов MPI;
•	MPI_MODE_NOCHECK — ни одна из ветвей не удерживает окно или не будет обращаться к окну во время удерживания блокировки. Это бывает полезно в тех случаях, когда взаимное исключение обеспечивается другими способами, но синхронизация все равно необходима.
Освобождение захваченного окна:
int MPI_Win_unlock(int rank, MPI_Win win);
Вызов этой функции завершает период RMA-доступа, начатый вызовом MPI_Win_lock (..., win). RMA-операции, вызванные во время этого периода, завершатся как в инициаторе, так и в адресате к моменту возврата из этого вызова. Блокировки используются, чтобы защитить обращения к заблокированному окну-адресату, на которое действуют RMA-вызовы, выданные между вызовами lock и unlock, и чтобы защитить локальные load/store обращения к заблокированному локальному окну, выполняемые между вызовами lock и unlock. Обращения, защищенные при помощи эксклюзивной блокировки, не будут пересекаться в пределах окна с другими обращениями к этому же окну, которое защищено блокировкой. Обращения, которые защищены совместной блокировкой, не будут пересекаться в пределах окна с обращениями, защищенными с помощью эксклюзивной блокировки, к одному и тому же окну.
6.9. Использование библиотеки МРЕ для анализа процессов взаимодействия ветвей программы
Библиотека МРЕ содержит программу визуализации лог-файлов Jumpshot-4 и набор функций, позволяющих создавать
125
эти лог-файлы и записывать в них «трассу» исполнения параллельной МР1-программы.
Приложение Jumpshot-4 позволяет представить работу параллельной программы в наиболее удобном для восприятия человека виде — визуализировать. Другими словами, программа Jumpshot-4 — это средство для визуализации трассы параллельной программы. Трасса — это журнал, содержащий информацию о ходе выполнения программы. Делает это она при помощи так называемых лог-файлов. Лог-файл — это файл трассы отмеченных во времени событий. Формат лог-файла, используемого Jumpshot-4, называется SLOG2.
Основной визуальный компонент Jumpshot-4 — канва с двумя измерениями: по горизонтали откладывается время, а по вертикали — процессы параллельной программы; канву можно масштабировать и прокручивать в горизонтальном и вертикальном направлениях.
Таким образом, каждая точка на канве может быть идентифицирована по двум числам — отметке времени и ID-линии времени.
Графические объекты, которые содержатся в SLOG-2 файле, рисуются на канве этой координатной системы. Эти объекты бывают двух типов: простые и составные. Простые объекты, как правило, являются основными объектами SLOG-2 файла. Есть деление этих объектов, основанное на их топологической структуре.
Формат SLOG-2 поддерживает три вида топологий: состояние, стрелка и событие. Состояние и стрелка определяются двумя точками на канве — это соответственно координаты: отметка времени и идентификатор (ID) линии времени. Начальный ID линии времени для состояния совпадает с конечным ID (т. е. этот объект относится к одному процессу), а для стрелки начальный и конечный ID могут различаться. Событие состоит только из одной точки на канве, т. е. имеет только одну линию времени и одно значение времени. Можно сказать, что состояние — это промежуток времени между двумя событиями.
Составные объекты имеют более сложную структуру и конструируются из набора простых объектов.
Лог-файлы и библиотека МРЕ
Исходные данные для программы Jumpshot содержатся в лог-файлах, которые можно получить путем включения в па-
126
раздельную программу вызовов функций библиотеки МРЕ (MultiprocessingEnvironment). Она содержит процедуры, которые облегчают написание, отладку и оценку эффективности MPI-программ. Рассмотрим структуру этой библиотеки чуть подробнее.
МРЕ-процедуры делятся на несколько категорий.
•	Параллельная графика (Parallel X graphics). Эти процедуры обеспечивают доступ всем процессам к разделяемому Х-дисплею. Они создают удобный вывод для параллельных программ, позволяют чертить текст, прямоугольники, круги, линии и т. д.
•	Регистрация (Logging). Библиотека МРЕ создает возможность легко получить лог-файл в каждом процессе и собрать их вместе по окончанию работы. Она также автоматически обрабатывает рассогласование и дрейф часов на множественных процессорах, если система не обеспечивает синхронизацию часов. Это библиотека для пользователя, который желает создать свои собственные события и программные состояния.
•	Последовательные секции (Sequential Sections). Иногда секция кода, которая выполняется на ряде процессов, должна быть выполнена только по одному процессу за раз в порядке номеров этих процессов. МРЕ имеет функции для такой работы.
•	Обработка ошибок (Error Handling). MPI имеет механизм, который позволяет пользователю управлять реакцией реализации на ошибки времени исполнения, включая возможность создать свой собственный обработчик ошибок.
Далее будут преимущественно рассмотрены средства регистрации (logging) для анализа эффективности параллельных программ. Анализ результатов регистрации производится после выполнения вычислений. Средства регистрации и анализа включают ряд профилирующих библиотек, утилитных программ и ряд графических средств.
Первая группа средств — профилирование. Библиотечные ключи обеспечивают собрание процедур, которые создают лог-файлы. Эти лог-файлы могут быть созданы вручную путем размещения в программе MPI обращений к МРЕ, или автоматически при установлении связи с соответствующими МРЕ- библиотеками (например, откомпилировав программу с библиотекой mpich2mpe.lib), или комбинацией этих двух методов.
127
В настоящее время МРЕ предлагает следующие три профилирующие библиотеки.
1.	Tracing Library (библиотека трассирования) — трассирует все MPI-вызовы. Каждый вызов предваряется строкой, которая содержит номер вызывающего процесса в MPI_COMM_WORLD и сопровождается другой строкой, указывающей, что вызов завершился. Большинство процедур SEND и RECIEVE также указывают значение count, tag и имена процессов, которые посылают или принимают данные.
2.	Animation Libraries (анимационные библиотеки) — простая форма программной анимации в реальном времени, которая требует процедур Х-окна.
3.	Logging Libraries (библиотеки регистрации) — самые полезные и широко используемые профилирующие библиотеки в МРЕ. Они формируют базис для генерации лог-файлов из пользовательских программ. Имеется три различных формата лог-файлов, допустимых в МРЕ. Формат CLOG содержит совокупность событий с единым отметчиком времени. Формат ALOG больше не развивается и поддерживается для обеспечения совместимости с ранними программами. И наиболее мощным является формат SLOG (для Scalable Logfile'), который может быть конвертирован из уже имеющегося CLOG-файла или получен прямо из выполняемой программы (для этого необходимо установить переменную среды MPE_LOG_FORMAT в SLOG).
Набор утилит в МРЕ включает конверторы лог-форматов (например, clog2slog), печать лог-файлов (slog_ print), оболочки средств визуализации лог-файлов, которые выбирают корректные графические средства для представления лог-файлов в соответствии с их расширениями.
Далее будут рассмотрены только библиотеки регистрации Logging Libraries. Оценки времени выполнения взаимодействий дают некоторое представление об эффективности программы. Но в большинстве случаев нужно подробно узнать, какова была последовательность событий, сколько времени было потрачено на каждую стадию вычисления и сколько времени занимает каждая отдельная операция передачи. Чтобы облегчить их восприятие, нужно представить их в графической форме.
Но для этого сначала нужно создать файлы событий со связанными временными отметками, затем исследовать их после окончания программы и только затем интерпретировать
128
их графически на рабочей станции. Такие файлы ранее уже названы лог-файлами. Способность автоматически генерировать лог-файлы является важной компонентой всех средств для анализа эффективности параллельных программ.
Опишем некоторые простые инструментальные средства для создания лог-файлов и их просмотра. Библиотека для создания лог-файлов отделена от библиотеки обмена сообщениями MPI. Просмотр лог-файлов независим от их создания, и поэтому могут использоваться различные инструментальные средства. Библиотека для создания лог-файлов МРЕ разработана таким образом, чтобы сосуществовать с любой МР1-реализацией и распространяется наряду с модельной версией MPI.
Чтобы создать файл регистрации, необходимо вызвать процедуру MPE_Log_event. Кроме того, каждый процесс должен вызвать процедуру MPE_Init_log, чтобы приготовиться к регистрации, и MPE_Finish_log, чтобы объединить файлы, сохраняемые локально при каждом процессе в единый лог-файл. МРЕ_ Stop_log используется, чтобы приостановить регистрацию, хотя таймер продолжает работать. MPE_Start_log возобновляет регистрацию.
Программист выбирает любые неотрицательные целые числа, желательные для типов событий; система не придает никаких частных значений типам событий. Чтобы присвоить событию свой номер, нужно вызвать функцию MPE_Log_get_ event_number. События рассматриваются, как не имеющие продолжительности. Чтобы измерить продолжительность состояния программы, необходимо, чтобы пара событий отметила начало и окончание состояния. Состояние определяется процедурой МРЕ-Describe _state, которая описывает начало и окончание типов событий. Процедура MPE_Describe_state также добавляет название состояния и его цвет на графическом представлении. Соответствующая процедура MPE_Describe_event обеспечивает описание события каждого типа.
Приведем прототипы описанных выше функций:
int MPE_Init_log (void)
int MPE_Start_log (void)
int MPE_Stop_log (void)
int MPE_Finish_log (char *logfilename)
int MPE Log get event number (void)
int MPE_Describe_state (int start, int end, char *name, char *color)
129
int MPE_Describe_event (int event, char *name)
int MPE_Log_event (int event, int intdata char *chardata)
Как уже было сказано, создать лог-файл можно и автоматически, откомпилировав проект с библиотекой mpich2mpe. lib (конечно, если при написании использовалась реализация MPICH2, аналогичная библиотека в MPICH, к сожалению, неизвестна). В этом случае в папке проекта создается файл *.clog2, который содержит всю трассу программы. В случае, если программист хочет какой-то «кусок» программы проанализировать отдельно или же просто создать свой лог-файл, имеет смысл использовать процедуры библиотеки МРЕ непосредственно в коде программы.
Рассмотрим, как это делается на примере простейшей программы (аналог «HelloWorld!»), которая печатает ранги всех процессов, на которых она выполнялась.
#include «stdafx.h»
#include "C:\Program Files\MPICH2\include\mpi.h"
#include "C:\Program Files\MPICH2\include\mpe.h"
int _tmain(int argc, char* argvf]) {
int ProcsNum,ProcRank;
int eventla,eventlb;
int event2a, event2b;
MPI_Init(&argc,&a rgv);
MPE_Init_log();
eventla = MPE Log get event numberQ;
eventlb = MPE Log get event numberQ;
event2a = MPE Log get event numberQ;
event2b = MPE Log get event numberQ;
MPE_Describe_state(eventla,eventlb,"GetSize","mage nta");
MPE_Describe_state(event2a,event2b,"GetRank","green");
MPE_Start_log();
MPE_Log_event(eventla,0,»start size»);
MPI_Comm_size(MPI_COMM_WORLD,&ProcsNum);
MPE_Log_event(eventlb,0,"finish size");
MPE_Log_event(event2a,0,"start rank");
MPI_Comm_rank(MPI_COMM_WORLD,&ProcRank);
MPE_Log_event(event2b,0,"finish rank");
printf("Hello from process %d\n",ProcRank);
MPE_Finish_log("logfilel");
MPI_Finalize(); return 0;
}
130
Вначале идет подготовка к регистрации (MPE_Init_log), после которой начинается блок описаний. В данной программе всего две функции MPI — получение размера системы и ранга процесса; предположим, мы хотим описать их обе. Для этого сначала определяются события. Поскольку каждое состояние (т. е. функция MPI) — это время между двумя событиями, то нам нужно определить четыре события. Для удобства лучше, чтобы названия событий, относящихся к одному состоянию, не сильно отличались (можно их назвать, например, eventla и eventlb). Функция MPE_Log_get_event_number присваивает каждому событию уникальный номер.
Далее описываются состояния: в данном случае состояние с названием «GetSize» расположено между событиями eventl а и eventlb и на канве в Jumpshot-4 будет отображаться малиновым цветом. Со вторым — аналогично.
Регистрация начинается путем вызова функции MPE_Start_ log. Считается, что события eventXa относятся к началу состояния, eventXb — к концу. Регистрируется определение размера системы (количества процессов). Для этого до и после соответствующей функции вызывается функция MPE_Log_event.
И, наконец, вызывается функция MPE_Finish_log, чтобы закончить регистрацию и закрыть лог-файл logfile.clog2.
Чтобы иметь возможность просмотреть созданный лог-файл в программе Jumpshot-4, нужно конвертировать его в logfile. slog2. Соответствующая утилита входит в состав приложения Jumpshot-4.
Контрольные вопросы и задания
1.	Какие виды виртуальных топологий можно реализовывать с использованием MPI?
2.	Что такое барьерная синхронизация в MPI? Какие еще виды синхронизации существуют в стандарте MPI-1?
3.	Перечислите основные группы функций стандарта MPI-1.
4.	Может ли MPI-программа, выполняющаяся на многоядерном узле, использовать все его ядра? И если да, то каким образом?
5.	Что такое пользовательская операция редукции?
6.	Чем различаются размер и протяженность типов данных MPI?
7.	Что такое буферизованная передача сообщений?
8.	Что такое коммуникатор? Какие проблемы решаются с использованием этого понятия?
131
9.	Сколько в MPI операций приема сообщений типа «точка-точка»? Перечислите и охарактеризуйте их.
10.	Что такое коллективное взаимодействие ветвей программы?
11.	Как путем вызова одной функции можно одновременно передать и принять блок данных?
12.	Перечислите предопределенные операции редукции данных MPI.
13.	Что такое виртуальная топология? Как она связана с физической топологией вычислительной сети/кластера/комплекса?
14.	Для чего в MPI-программах можно использовать производные типы данных?
15.	Каково назначение библиотеки МРЕ?
16.	Что понимается под операциями сдвига и циклического сдвига данных? Какие функции обычно используются для их реализации?
17.	Что такое исключающая редукция (MPI_Exscan), как она выполняется?
18.	Может ли ветвь параллельной программы с помощью MPI передавать сообщения сама себе?
19.	Какие операции рассылки и сбора данных реализованы в MPI?
20.	Должны ли абсолютно все ветви параллельной программы принимать участие в коллективном взаимодействии?
21.	Сколько конструкторов производных типов реализовано в MPI? Перечислите и охарактеризуйте их.
22.	Чем различаются топологии «тор» и «решетка»?
23.	Какие способы синхронизации используются при удаленном доступе к памяти?
24.	Что такое «окно» при удаленном доступе памяти? Как создать окно?
25.	Что такое синхронизация с блокировками при удаленном доступе к памяти? Какие еще существуют виды синхронизации?
26.	Какие функции могут быть использованы для модификации состояния памяти «окна» другой ветви?
27.	Перечислите основные группы функций библиотеки МРЕ.
28.	Что такое заборная синхронизация?
29.	С помощью каких функций выполняется создание новых групп ветвей?
30.	Что такое синхронизация с пассивным адресатом?
31.	С помощью каких функций выполняются захват/освобождение окна при удаленном доступе к памяти?
Рекомендуемая литература
1.	Воеводин, В. В. Параллельные вычисления / В. В. Воеводин, Вл. В. Воеводин. — Санкт-Петербург : БХВ-Петербург, 2004.
2.	Воеводин В. В. Вычислительная математика и структура алгоритмов / В. В. Воеводин — Москва : Изд-во МГУ, 2006.
3.	Таненбаум Э. Архитектура компьютера / Э. Таненбаум. — 5-е изд. — Санкт-Петербург : Питер, 2007. 848 с.
4.	Новожилов, О. П. Архитектура ЭВМ и систем. В 2 ч. Часть 2 : учебное пособие для академического бакалавриата / О. П. Новожилов. — Москва : Издательство Юрайт, 2019.
5.	Баденко В. Л. Высокопроизводительные вычисления : учебное пособие. — Санкт-Петербург : Изд-во Политехи, ун-та, 2010. — 180 с.
6.	Шпаковский Г. И. Реализация параллельных вычислений: кластеры, многоядерные процессоры, грид, квантовыекомпью-теры. — Минск, БГУ, 2010 г., 155 с.
7.	Рыбалъченко, М. В. Архитектура информационных систем : учебное пособие для вузов / М. В. Рыбальченко. — Москва : Издательство Юрайт, 2019.
8.	Бабичев, С. Л. Распределенные системы : учебное пособие для вузов / С. Л. Бабичев, К. А. Коньков. — Москва : Издательство Юрайт, 2019.
9.	Миллер, Р. Последовательные и параллельные алгоритмы. Общий подход / Р. Миллер, Л. Боксер. — Москва : Бином. Лаборатория знаний, 2015.
10.	Лупин, С. А. Технологии параллельного программирования / С. А. Лупин, М. А. Посыпкин. — Москва : ИД Форум, 2011.
11.	Антонов, А. С. Параллельное программирование с использованием технологии ОрепМР / А. С. Антонов. — Москва : Изд-во МГУ, 2009.
12.	Левин, М. П. Параллельное программирование с использованием ОрепМР / М. П. Левин. — Москва : Бином. Лаборатория знаний, 2016.
133
13.	Официальный сайт ОрепМР Architecture Review Board. — URL: http://openmp.org/ (дата обращения 03.12.2019).
14.	Боресков, А. В. Основы CUDA / А. В. Боресков. — URL: http: //www. steps3d .narod. ru/tutorials/cuda-tutorial.html (дата обращения 03.12.2019).
15.	Сандерс Д., Кондрат Э. Технология CUDA в примерах. Введение в программирование графических процессоров. — Москва : ДМК Пресс, 2011 г.
16.	Боресков, А. В. Параллельные вычисления на GPU. Архитектура и программная модель CUDA / А. В. Боресков. — Москва : Изд-во МГУ, 2015.
17.	Борисов, Е. С. Технология параллельного программирования CUDA / Е. С. Борисов. — URL: http://mechanoid.kiev.ua/ parallel-cuda.html (дата обращения 06.12.2019).
18.	Официальный сайт разработчиков NVIDIA CUDA. — URL: https://developer.nvidia.com/cuda-zone (дата обращения 06.12.2019).
19.	Антонюк В. A. OpenCL. Открытый язык для параллельных программ / В. А. Антонюк. — Москва : Изд-во МГУ, 2017
20.	Казённое А. М. Основы технологии CUDA и OpenCL / А. М. Казеинов. — Москва : Изд-во МФТИ, 2013.
21.	Копысов С. П., Новиков А. К. Промежуточное программное обеспечение параллельных вычислений : учебное пособие / С. П. Копысов, А. К. Новиков. — Ижевск : Изд-во Удмуртского государственного университета, 2012.
22.	Официальный сайт Khronos group. OpenCL. — URL: https://www.khronos.org/opencl/ (дата обращения 06.12.2019).
23.	Корнеев, В. Д. Параллельное программирование в MPI / В. Д. Корнеев. — Москва ; Ижевск : Институт компьютерных исследований, 2003.
24.	Антонов, А. С. Параллельное программирование с использованием технологии MPI / А. С. Антонов. — Москва : Изд-во МГУ, 2004.
25.	Сайт форума разработчиков MPI. — URL: https://www. mpi-forum.org/ (дата обращения 09.12.2019).
26.	Официальный сайт разработчиков MPICH. — URL: http://www. mcs.anl.gov/project/mpich-high-performance-portable-implementation-mpi/ (дата обращения 03.12.2019).
27.	Богачев, К Ю. Основы параллельного программирования / К. Ю. Богачев — Москва : Бином. Лаборатория знаний, 2015.
134
28.	Гергелъ, В. П. Высокопроизводительные вычисления для многопроцессорных многоядерных систем / В. П. Гергель. — Москва : Изд-во МГУ, 2011.
29.	Тормасов, А. Г. Параллельное программирование многопоточных систем с разделяемой памятью / А. Г. Тормасов. — Москва : Изд-во МФТИ, 2014.