Введение. Вычислительные системы. Идеология OpenMP. Параллельная обработка в конструкциях, отличных от циклов

OpenMP позволяет легко создавать параллельные программы для систем с общей памятью (много­ядерных и многопроцессорных). В этой статье я расскажу о том, как включить OpenMP в самой рас­пространённой среде программирования - Visual Studio. Согласно официальной версии Micro­soft , OpenMP поддерживается только в версиях Professional среды разработки Visual Studio 2005/2008/2010. Однако бесплатная Visual Studio Express обладает тем же компилятором, что и версия Professional. Поэтому после небольшой «доработки напильником» параллельные OpenMP-программы будут компилироваться и, главное, работать даже в Visual Studio Express.

OpenMP от Microsoft реализован посредством следующих компонентов:

1. компилятор C++, входящий в состав Visual Studio;
2. заголовочный файл omp.h;
3. библиотеки стадии компиляции: vcomp.lib и vcompd.lib (последняя используется для отлад­ки);
4. библиотеки времени выполнения: vcomp90.dll и vcomp90d.dll. Цифры в названии могут разли­чаться: в Visual Studio 2005 вместо 90 стоят цифры 80.

В бесплатной Visual Studio Express перечисленные библиотеки отсутствуют.

OpenMP и Visual Studio Express

Если вы хотите создавать параллельные OpenMP-программы под Windows, то самый удобный спо­соб - это воспользоваться Visual Studio 2005/2008/2010 Professional. Напоминаю, что она бес­платна для студентов и аспирантов . Кроме того, вы можете купить Visual Studio Professional за $600 (конечно, существует ещё и третий вариант, но мы о нём не будем говорить).

Если на вашем компьютере установлена версия Professional - переходите к следующему разделу статьи. В этом разделе рассмотрим случай, когда по каким-либо причинам вы вынуждены использо­вать Visual Studio Express.

Visual Studio Express - это бесплатная урезанная версия Visual Studio. Нас будет интересовать вер­сия 2008 года. Скачать её можно отсюда: http://www.microsoft.com/exPress/download/ . Естествен­но, программировать будем на C++, поэтому выби­райте Visual C++ 2008.

Программа установки будет загружать данные из Интернета (примерно 100 мегабайт), поэтому вы можете сэкономить немного трафика, отключив установку Microsoft Silverlight и Microsoft SQL Server, если они вам не требуются.

После установки Visual Studio Express нам понадо­бится добавить в неё OpenMP-компоненты. Легаль­ный бесплатный способ сделать это - установить Windows SDK for Windows Server 2008 and .NET Framework 3.5 . Во время установки этого пакета программ будет произведено обновление Visual Studio. Процесс обновления не смотрит, какая именно версия Visual Studio у вас установлена (Express или Professional), поэтому во время установ­ки будут «случайно» добавлены недостающие компоненты.

Так как мы ставим Windows SDK только ради OpenMP, то нам не нужен тот гигабайт документации, который идёт в комплекте. Рекомендую оставить только следующие элементы:

Рисунок 1. Необходимые нам компоненты SDK

К сожалению, в состав SDK не входит библиотека vcomp90d.dll, поэтому на данный момент в Visual Studio Express вы сможете запускать только OpenMP-программы, откомпилированные в ре­жиме Release. Я нашёл способ обойти и это ограничение, об этом читайте далее (раздел «Отладка OpenMP-программы в Visual Studio Express»).

Использование OpenMP в Visual Studio

После того, как вы выполнили шаги, описанные в предыдущем разделе, уже не важно, какой верси­ей Visual Studio вы пользуетесь. Покажу шаг за шагом, как создать проект с поддержкой OpenMP в этой среде разработки. Прежде всего, нужно запустить Visual Studio, и выбрать File→New→Project... Появится окно создания проекта. Выберите тип проекта «Win32», шаблон - «Win32 Console Application». Введите осмысленное имя проекта, выберите папку для хранения проекта, уберите галочку «Create directory for solution»:

Рисунок 2. Окно создания проекта

Нажмите кнопку «OK», появится окно настройки будущего проекта. Выберите вкладку «Application Settings», и включите галку «Empty project»:

Рисунок 3. Окно настройки будущего проекта

По нажатию кнопки «Finish» проект будет создан. Никаких видимых изменений в главном окне Visual Studio не произойдёт. Только имя проекта в заголовке окна как бы говорит нам о том, что мы ра­ботаем с проектом.

Теперь нажмите Project→Add New Item, появится окно добавления элементов в проект. Добавьте.cpp-файл в проект:

Рисунок 4. Окно добавления элементов в проект

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

#include #include using namespace std; int main(int argc, char **argv) { int test(999); omp_set_num_threads(2); #pragma omp parallel reduction(+:test) { #pragma omp critical cout << "test = " << test << endl; } return EXIT_SUCCESS; } Листинг 1. Простейшая программа, использующая OpenMP

Запустите программу, нажав Debug→Start Without Debugging. Если всё было сделано правильно, программа откомпилируется (если спросит вас, компилировать ли, нажмите «Yes»), затем запу­стится и выведет test = 999:

Рисунок 5. Результат работы программы из листинга 1

«Как же так?! - скажете вы - Ведь программа должна была вывести ноль, причём дважды!». Дело в том, что OpenMP ещё не включен, и поэтому соответствующие директивы были проигнорированы компилятором.

Для включения OpenMP нажмите Project→OMP Properties (OMP - имя проекта из моих примеров). Слева вверху появившегося окна выберите «All Configurations» и в разделе Configuration Properties→C/C++→Language включите «OpenMP Support»:

Рисунок 6. Включаем OpenMP в свойствах проекта

После этого снова запустите программу, нажав Debug→Start Without Debugging. На этот раз про­грамма выведет test = 0 дважды:

Рисунок 7. Результат работы программы из листинга 1 с включённым OpenMP

Ура! OpenMP работает.

Примечание. Если вы используете Visual Studio Express, то выберите текущую конфигурацию «Release», иначе работать не будет (читайте далее):

Рисунок 8. Выбор текущей конфигурации

Как было сказано ранее, даже после установки Windows SDK у нас не будет в наличии необходимой для отладки библиотеки vcomp90d.dll, поэтому мы пока не можем отлаживать OpenMP програм­му в Visual Studio Express. Простое копирование имеющейся библиотеки vcomp90.dll и переимено­вание её в vcomp90d.dll не сработает, ибо не совпадёт контрольная сумма и версия, указанные во встраиваемом в exe-файл манифесте. Поэтому будем «копать» с противоположной стороны.

При компиляции в конфигурации «Debug» («Отладка»), заголовочный файл omp.h требует библио­теку vcompd.lib (она у нас имеется), которая, в свою очередь, требует vcomp90d.dll (отсутству­ет). Лицензия не позволяет нам использовать в приложениях модифицированные заголовочные файлы от Microsoft, поэтому вместо модификации omp.h включим его в нашу программу следую­щим образом, чтобы он не догадался о включённом режиме отладки:

#include #ifdef _DEBUG #undef _DEBUG #include #define _DEBUG #else #include #endif using namespace std; int main(int argc, char **argv) { int test(999); omp_set_num_threads(2); #pragma omp parallel reduction(+:test) { #pragma omp critical cout << "test = " << test << endl; } return EXIT_SUCCESS; } Листинг 2. Включаем omp.h «хитрым» способом

Приведённого действия не достаточно для того, чтобы всё работало (пока мы исправили лишь ма­нифест, встраиваемый в программу). Дело в том, что Visual Studio в режиме отладки по прежнему автоматически (из-за включённого OpenMP) прилинковывает vcompd.lib, требующую vcomp90d.dll. Чтобы это исправить, снова зайдите в настройки проекта (Project→OMP Properties), выберите на этот раз Configuration: «Debug». В разделе Configuration Properties→Linker→Input укажите, что vcompd.lib прилинковывать не нужно, а vcompd.lib - нужно:

Рисунок 9. Заменяем библиотеку в свойствах проекта

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

Рисунок 10. Точка останова

После этого запустите программу в режиме отладки: Debug→Start Debugging (не забудьте вер­нуть текущую конфигурацию «Debug», см. рисунок 8). Программа запустится - и сразу же остано­вится на точке останова. Во вкладке «Threads» мы видим, что программа действительно работает, используя два потока:

Рисунок 11. Отладка OpenMP-программы в Visual Studio Express

Библиотека активно развивается, в настоящий момент актуальный стандарт версии 4.5 (выпущен в 2015 году), однако уже идет разработка версии 5.0. В тоже время, компилятор Microsoft C++ поддерживает только версию 2.0, а gcc — версию 4.0.

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

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

• в книгах стараются описать все возможности OpenMP, значительная часть которых является атавизмом. В старых версия библиотеки это было нужно, но сейчас использование этих возможностей делает ваш код опаснее и повышает порог вхождения (другой программист потратит много времени чтобы разобраться в нем). Я кратко описал некоторые из таких возможностей в заметке « » и в «учебнике» про них ничего нет. За счет этого мой учебник получился короче и не загружает читателя устаревшей информацией;
• в книгах мало «реальных» примеров — поиск минимума/максимума не в счет. Я предлагаю своим читателям обратиться за такими примерами на . На форуме можно взять исходный код программ, которые отлажены и работают, при этом к каждой программе прилагается подробный разбор решения и подхода к распараллеливанию. Иногда можно найти несколько реализаций. Примеры будут дополняться.

Вычислительные системы. Идеология OpenMP

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

Для использования библиотеки OpenMP вам необходимо подключить заголовочный файл "omp.h" , в а также добавить опцию сборки -fopenmp (для компилятора gcc) или установить соответствующий флажок в настройках проекта (для Visual Studio). После запуска программы создается единственный процесс, который начинается выполняться как и обычная последовательная программа. Встретив параллельную область (задаваемую директивой #pragma omp parallel) процесс порождает ряд потоков (их число можно задать явно, однако по умолчанию будет создано столько потоков, сколько в вашей системе вычислительных ядер). Границы параллельной области выделяются фигурными скобками, в конце области потоки уничтожаются. Схематично этот процесс изображен на рисунке 2:

рис.2 директива omp parallel

Черными линиями на рисунке показано время жизни потоков, а красными — их порождение. Видно, что все потоки создаются одним (главным) потоком, который существует все время работы процесса. Такой поток в OpenMP называется master , все остальные потоки многократно создаются и уничтожаются. Стоит отметить, что директивы parallel могут быть вложенными, при этом в зависимости от настроек (изменяются функцией omp_set_nested) могут создаваться вложенные потоки.

OpenMP может использоваться на кластерной архитектуре, но не самостоятельно, т.к. загрузить весь кластер она не может (для этого нужно создавать процессы, использовать инструменты типа MPI ). Однако, если узел кластера многоядерный — то использование OpenMP может существенно поднять эффективность. Такое приложение будет являться «гибридным», в этой статье я не буду вдаваться в тему, тем более что по ней написано много хороших материалов .

Синхронизация — критические секции, atomic, barrier

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

#include "omp.h" #include int main() { int value = 123; #pragma omp parallel { value++; #pragma omp critical { std::cout << value++ << std::endl; } } }

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

Проблема гонки потоков OpenMP

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

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

Для ряда операций более эффективно использовать директиву atomic , чем критическую секцию. Она ведет себя также, но работает чуть быстрее. Применять ее можно для операций префиксного/постфиксного инкремента/декремента и операции типа X BINOP = EXPR , где BINOP представляет собой не перегруженный оператор +, *, -, /, &, ^, |, <<, >> . Пример использования такой директивы:

#include "omp.h" #include int main() { int value = 123; #pragma omp parallel { #pragma omp atomic value++; #pragma omp critical (cout) { std::cout << value << std::endl; } } }

Тут же показана возможность именования критических секций — очень рекомендую ей пользоваться всегда. Дело в том, что все безымянные секции рассматриваются как одна (очень большая) и если вы не дадите им явно имена, то только в одной из этих секций в один момент времени будет один поток. Остальные будут ждать. Имя секции должно подсказывать программисту к какому виду ресурса она относится — в приведенном примере таким ресурсом является поток вывода на экран (cout).

Несмотря на то, что каждая операция над общими данными в последнем примере размещена в критической секции или является атомарной, в нем есть проблема, т.к. порядок выполнения этих операций по прежнему не определен. Запустив программу раз 20 мне удалось получить на экране не только "125 125" , но и "124 125" . Если мы хотим чтобы сначала каждый поток увеличил значение, а затем они вывели их на экран — можем использовать директиву barrier:

#pragma omp parallel { #pragma omp atomic value++; #pragma omp barrier #pragma omp critical (cout) { std::cout << value << std::endl; } }

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

Разделение задач между потоками

Параллельный цикл

Самый популярный способ распределения задач в OpenMP — параллельный цикл. Не секрет, что программы почти всю свою жизнь проводят выполняя циклы, при этом если между итерациями цикла нет зависимостей — то цикл называется векторизуемым (его итерации можно поделить между потоками и выполнить независимо друг от друга). В статье « » я поверхностно описал эту конструкцию на примерах вычисления суммы элементов массива и численного интегрирования, не буду повторяться — рекомендую вам пройти по ссылке, т.к. дальше описаны более «интересные» аспекты параллельного цикла.

Параллельный цикл позволяет задать опцию schedule , изменяющую алгоритм распределения итераций между потоками. Всего поддерживается 3 таких алгоритма. Далее полагаем, что у нас p потоков выполняют n итераций:

Опции планирования:

• schedule(static) — статическое планирование. При использовании такой опции итерации цикла будут поровну (приблизительно) поделены между потоками. Нулевой поток получит первые \(\frac{n}{p}\) итераций, первый — вторые и т.д.;
• schedule(static, 10) — блочно-циклическое распределение итераций. Каждый поток получает заданное число итераций в начале цикла, затем (если остались итерации) процедура распределения продолжается. Планирование выполняется один раз, при этом каждый поток «узнает» итерации которые должен выполнить;
• schedule(dinamic), schedule(dynamic, 10) — динамическое планирование. По умолчанию параметр опции равен 1. Каждый поток получает заданное число итераций, выполняет их и запрашивает новую порцию. В отличии от статического планирования, выполняется многократно (во время выполнения программы). Конкретное распределение итераций между потоками зависит от темпов работы потоков и трудоемкости итераций;
• schedule(guided), schedule(guided, 10) — разновидность динамического планирования с изменяемым при каждом последующем распределении числе итераций. Распределение начинается с некоторого начального размера, зависящего от реализации библиотеки до значения, задаваемого в опции (по умолчанию 1). Размер выделяемой порции зависит от количества еще нераспределенных итераций

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

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

Int sum_arr(int *a, const int n) { int sum = 0; #pragma omp parallel reduction (+: sum) { #pragma omp for for (int i = 0; i < n; ++i) sum += a[i]; } return sum; }

Выглядит это красиво, но на самом деле в каждом потоке создается локальная переменная для хранения суммы части массива (вычисление которой назначено текущему потоку), ей присваивается значение 0 (т.к. редукция с оператором +). Каждый поток вычисляет сумму, но необходимо ведь сложить все эти значение чтобы получить окончательный результат? — Делается это с помощью критической секции или атомарной операции примерно следующим образом:

Int sum_arr(int *a, const int n) { int sum = 0; #pragma omp parallel { int local_sum = 0; #pragma omp for for (int i = 0; i < n; ++i) local_sum += a[i]; #pragma omp atomic sum += local_sum; } return sum; }

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

Параллельные задачи (parallel tasks)

Параллельные задачи — это более гибкий механизм, чем параллельный цикл. Параллельный цикл описывается внутри параллельной области, при этом могут возникнуть проблемы. Например, мы написали параллельную функцию вычисления суммы элементов одномерного массива, и нашу функцию решили применить для вычисления суммы элементов матрицы, но сделать это также параллельно. Получится вложенный параллелизм. Если (теоретически) наш код запущен на 8 ядрах — то фактически будет создано 64 потока. Ну а если кому-нибудь придет в голову идея делать параллельно еще что-нибудь?

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

Проблема параллельного цикла в том, что число создаваемых потоков зависит от того какие функции распараллелены и как они друга друга вызывают. Очень сложно все это отслеживать и, тем более, поддерживать. Решение проблемы — параллельные задачи, которые не создают поток, а лишь выполняют добавление задачи в очередь, освободившийся поток выбирает задачу из пула. Я описывал этот механизм в статье
« » и не буду повторяться (рекомендую прочитать материал по ссылке — в статье рассматривается возможность распараллеливания рекурсивных функций с помощью механизма задач). Отмечу лишь то, что в параллельные задачи были предложены в стандарте OpenMP 3.0 (в 2008 году) поэтому их поддержка отсутствует в Microsoft C++. Кроме того, в свежем стандарте OpenMP 4.5 была предложена конструкция taskloop , за счет которой использовать параллельные задачи для распараллеливания циклов теперь также удобно как и параллельный цикл.

Параллельные секции

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

For (int i = 1; i < n; ++i) a[i] = a+1;

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

#pragma omp parallel { #pragma omp sections { #pragma omp section { for (int i = 1; i < n; ++i) a[i] = a+1; } #pragma omp section { for (int i = 1; i < n; ++i) b[i] = b+1; } } }

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

Заключение и дополнительная литература

Тут учебнику конец… Если хотите посмотреть больше примеров исходников с OpenMP или у вас есть вопросы по решению конкретных задач — загляните на наш форум. Если же у вас возникли вопросы по тексту статьи — пишите в комментариях.

1. Поддержка OpenMP в Lazarus/Freepascal: http://wiki.lazarus.freepascal.org/OpenMP_support
2. OpenMP для Java: http://www.omp4j.org/
3. Антонов А.С. Технологии параллельного программирования MPI и OpenMP: Учеб. пособие. Предисл.: В.А.Садовничий. – М.: Издательство Московского университета, 2012.-344 с.-(Серия “Суперкомпьютерное образование”). ISBN 978-5-211-06343-3.
4. Гергель В.П. Высокопроизводительные вычисления дл многоядерных многопроцессорных систем. Учебное пособие – Нижний Новгород; Изд-во ННГУ им. Н.И.Лобачевского, 2010
5. : https://сайт/archives/1150
6. Параллельные задачи (tasks) OpenMP:
7. : https://сайт/forums/forum/programming/parallel_programming/openmp_library
8. : href=»https://сайт/forums/topic/openmp-problems
9. : https://сайт/forums/topic/matrix-triangulation_cplusplus
10. Профилировка гибридных кластерных приложений MPI+OpenMP: https://habrahabr.ru/company/intel/blog/266409/
11. : https://сайт/forums/topic/openmp-cramer-method_cplusplus
12. 32 подводных камня OpenMP при программировании на Си++: https://www.viva64.com/ru/a/0054/
13. OpenMP и статический анализ кода: https://www.viva64.com/ru/a/0055/

Инфраструктура OpenMP позволяет эффективно реализовывать технологии параллельного программирования на C, C++ и Fortran. Компилятор GNU Compiler Collection (GCC) версии 4.2 поддерживает спецификацию OpenMP 2.5, а GCC версии 4.4 – самую последнюю спецификацию OpenMP 3. Другие компиляторы, включая Microsoft® Visual Studio, также поддерживают OpenMP. Эта статья научит вас использовать прагмы компилятора OpenMP; также она содержит информацию о некоторых API-интерфейсах OpenMP и раскрывает некоторые приемы параллельных вычислений с использованием OpenMP. Во всех примерах этой статьи используется компилятор GCC 4.2.

Начало работы

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

Создаем первое OpenMP-приложение

Начнем с написания простого приложения Hello, World! , содержащего дополнительную прагму. Код этого приложения представлен в листинге 1.

Листинг 1. Программа "Hello World", написанная с использованием OpenMP

#include int main() { #pragma omp parallel { std::cout << "Hello World!\n"; } }

После компиляции и запуска этого кода с помощью g++ вы увидите на экране надпись Hello, World! . Теперь перекомпилируем код с опцией -fopenmp . Результат работы программы представлен в листинге 2.

Листинг 2. Компиляция и запуск кода с использованием опции -fopenmp

tintin$ g++ test1.cpp -fopenmp tintin$ ./a.out Hello World! Hello World! Hello World! Hello World! Hello World! Hello World! Hello World! Hello World!

Что же произошло? Когда вы используете опцию компилятора -fopenmp , в дело вступает директива #pragma omp parallel . В процессе компиляции внутренние механизмы GCC создают столько параллельных потоков, сколько могут выполняться при условии оптимальной загрузки системы (в зависимости от конфигураций аппаратного обеспечения и операционной системы), при этом в каждом создаваемом потоке выполняется код, заключенный в блоке после прагмы. Такое поведение называется неявным распараллеливанием , а ядро OpenMP состоит из набора мощных прагм, избавляющих вас от написания множества типовых фрагментов кода (для интереса вы можете сравнить приведенный кода с реализацией этой же задачи с помощью POSIX-потоков ). Я использую компьютер с процессором Intel® Core i7 с 4 физическими ядрами по 2 логических ядра в каждом, что объясняет результаты, полученные в листинге 2 (8 потоков = 8 логических ядер).

Возможности OpenMP parallel

Количеством потоков можно легко управлять с помощью прагмы с аргументом num_threads . Ниже представлен код из листинга 1 с заданным количеством потоков (5 потоков):

Листинг 3. Управление количеством потоков с помощью num_threads

#include int main() { #pragma omp parallel num_threads(5) { std::cout << "Hello World!\n"; } }

Вместо аргумента num_threads можно воспользоваться альтернативным способом для задания количества потоков выполнения кода. Здесь мы подошли к рассмотрению первого API-интерфейса OpenMP под названием omp_set_num_threads . Эта функция определяется в файле заголовка omp.h. Для выполнения кода из листинга 4 не требуется использовать какие-либо дополнительные библиотеки, – просто используйте опцию -fopenmp .

Листинг 4. Точное управление потоками с помощью omp_set_num_threads

#include #include int main() { omp_set_num_threads(5); #pragma omp parallel { std::cout << "Hello World!\n"; } }

Наконец, для управления работой OpenMP можно использовать внешние переменные среды. Можно исправить код из листинга 2 и просто напечатать фразу Hello World! шесть раз, задав для переменной OMP_NUM_THREADS значение 6, как показано в листинге 5.

Листинг 5. Переменные среды для настройки OpenMP

tintin$ export OMP_NUM_THREADS=6 tintin$ ./a.out Hello World! Hello World! Hello World! Hello World! Hello World! Hello World!

Мы рассмотрели три стороны OpenMP: прагмы компилятора, API-интерфейсы времени выполнения и переменные среды. Что произойдет, если использовать переменные среды вместе с API-интерфейсами времени выполнения? API-интерфейсы имеют более высокий приоритет.

Практический пример

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

Листинг 6. Последовательная обработка в цикле for

int main() { int a, b; // ... код инициализации массивов a и b; int c; for (int i = 0; i < 1000000; ++i) c[i] = a[i] * b[i] + a * b; // ... выполняем некоторые действия с массивом c }

Очевидно, что цикл for можно распараллелить и обрабатывать сразу несколькими ядрами процессора, поскольку вычисление значения любого элемента c[k] никак не зависит от остальных элементов массива c . В листинге 7 показано, как можно это сделать с помощью OpenMP.

Листинг 7. Параллельная обработка в цикле for с помощью прагмы parallel for

int main() { int a, b; // ... код для инициализации массивов a и b; int c; #pragma omp parallel for for (int i = 0; i < 1000000; ++i) c[i] = a[i] * b[i] + a * b; // ... выполняем некоторые действия с массивом c }

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

Листинг 8. Пример с использованием API-функции omp_get_wtime

#include #include #include #include int main(int argc, char *argv) { int i, nthreads; clock_t clock_timer; double wall_timer; double c; for (nthreads = 1; nthreads <=8; ++nthreads) { clock_timer = clock(); wall_timer = omp_get_wtime(); #pragma omp parallel for private(i) num_threads(nthreads) for (i = 0; i < 1000000; i++) c[i] = sqrt(i * 4 + i * 2 + i); std::cout << "threads: " << nthreads << " time on clock(): " << (double) (clock() - clock_timer) / CLOCKS_PER_SEC << " time on wall: " << omp_get_wtime() - wall_timer << "\n"; } }

В листинге 8 мы измеряем время выполнения внутреннего цикла for , увеличивая при этом количество потоков. API-функция omp_get_wtime возвращает затраченное фактическое время (в секундах), прошедшее с начала заданной точки отсчета. Таким образом, значение omp_get_wtime() - wall_timer возвращает фактическое время выполнения цикла for . Системный вызов clock() используется для оценки времени, затраченного центральным процессором на выполнение всей программы, т. е. прежде чем получить итоговый результат, мы суммируем все эти временные интервалы с учетом потоков. На моем компьютере с процессором Intel Core i7 я получил результаты, приведенные в листинге 9.

Листинг 9. Статистика выполнения внутреннего цикла for

threads: 1 time on clock(): 0.015229 time on wall: 0.0152249 threads: 2 time on clock(): 0.014221 time on wall: 0.00618792 threads: 3 time on clock(): 0.014541 time on wall: 0.00444412 threads: 4 time on clock(): 0.014666 time on wall: 0.00440478 threads: 5 time on clock(): 0.01594 time on wall: 0.00359988 threads: 6 time on clock(): 0.015069 time on wall: 0.00303698 threads: 7 time on clock(): 0.016365 time on wall: 0.00258303 threads: 8 time on clock(): 0.01678 time on wall: 0.00237703

Хотя время центрального процессора (time on clock) во всех случаях получилось примерно одинаковым (как и должно быть, не принимая во внимание дополнительное время, затраченное на создание потоков и контекстного переключателя), интересующее нас фактическое время (time on wall) постоянно уменьшалось при увеличении количества потоков, которые предположительно выполнялись параллельно отдельными процессорными ядрами. Итак, сделаем последнее замечание относительно синтаксиса прагмы: #pragma parallel for private(i) означает, что переменная цикла i рассматривается как локальная память потока; каждый поток содержит свою копию этой переменной. Локальная переменная потока не инициализируется.

Критические участки кода в OpenMP

Разумеется, вы понимаете, что нельзя полностью доверить OpenMP автоматически обрабатывать критические участки кода, не так ли? Конечно, вам не придется явным образом создавать взаимные исключения (мьютексы), тем не менее критические участки необходимо указывать. Синтаксис приводится в следующем примере:

#pragma omp critical (optional section name) { // 2 потока не могут выполнять этот блок кода одновременно }

Код, следующий после директивы pragma omp critical , может выполняться в заданный момент времени только в одном потоке. Кроме того, имя раздела optional section name является глобальным идентификатором, и критические участки с одинаковым идентификатором не могут обрабатываться одновременно двумя потоками. Рассмотрим код в листинге 10.

Листинг 10. Несколько критических участков с одинаковыми именами

#pragma omp critical (section1) { myhashtable.insert("key1", "value1"); } // ... здесь содержится какой-то другой код #pragma omp critical (section1) { myhashtable.insert("key2", "value2"); }

Посмотрев на код в этом листинге, можно предположить, что две операции вставки данных в хэш-таблицу никогда не будут выполняться одновременно, поскольку имена критических участков совпадают. Это немного отличается от того, к чему вы привыкли, обрабатывая критические участки в pthreads, для которых характерно использование большого количества блокировок (что может приводить к лишним сложностям).

Блокировки и мьютексы в OpenMP

Интересно, что OpenMP содержит свои версии мьютексов (в конце концов, OpenMP – это не только прагмы). Итак, познакомьтесь с типом omp_lock_t , определенным в заголовочном файле omp.h. Обычные мьютексные операции в стиле pthread-потоков содержат значение true, даже если имена API-функций одинаковы. Вот пять API-функций, о которых необходимо знать:

• omp_init_lock : эта API-функция должна использоваться первой при обращении к типу omp_lock_t и предназначена для инициализации. Необходимо отметить, что сразу после инициализации блокировка будет находиться в исходном (unset) состоянии.
• omp_destroy_lock : уничтожает блокировку. В момент вызова этой API-функции блокировка должна находиться в исходном состоянии; это означает, что нельзя вызвать функцию omp_set_lock , а затем уничтожить блокировку.
• omp_set_lock : устанавливает omp_lock_t , т. е. активирует мьютекс. Если поток не может установить блокировку, он продолжает ожидать до тех пор, пока такая возможность не появится.
• omp_test_lock : пытается установить блокировку, если она доступна; возвращает 1 в случае успеха и 0 в случае неудачи. Это функция является неблокирующей , т. е. она не заставляет поток ожидать установления блокировки.
• omp_unset_lock : сбрасывает блокировку в исходное состояние.

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

Листинг 11. Усовершенствование однопоточной очереди с помощью OpenMP

#include #include "myqueue.h" class omp_q: public myqueue { public: typedef myqueue base; omp_q() { omp_init_lock(&lock); } ~omp_q() { omp_destroy_lock(&lock); } bool push(const int& value) { omp_set_lock(&lock); bool result = this->base::push(value); omp_unset_lock(&lock); return result; } bool trypush(const int& value) { bool result = omp_test_lock(&lock); if (result) { result = result && this->base::push(value); omp_unset_lock(&lock); } return result; } // likewise for pop private: omp_lock_t lock; };

Вложенные блокировки

Другими типами блокировок в OpenMP являются различные варианты блокировки omp_nest_lock_t . Они похожи на omp_lock_t , но имеют дополнительное преимущество: такая блокировка может быть включена удерживающим ее потоком несколько раз. Каждый раз, когда вложенная блокировка активируется удерживающим ее потоком с помощью omp_set_nest_lock , увеличивается внутренний счетчик блокировок. Блокировка освобождается удерживающим потоком тогда, когда один или несколько вызовов omp_unset_nest_lock снижают значение внутреннего счетчика блокировок до 0 . Для работы с omp_nest_lock_t используются следующие API-функции:

• omp_init_nest_lock(omp_nest_lock_t*) : устанавливает внутренний счетчик вложенности в 0 .
• omp_destroy_nest_lock(omp_nest_lock_t*) : уничтожает блокировку. Вызов этой API-функции для блокировки со значением счетчика, отличного от нуля, приводит к непредсказуемым результатам.
• omp_set_nest_lock(omp_nest_lock_t*) : аналогична функции omp_set_lock за исключением того, что удерживающий поток может вызывать ее несколько раз.
• omp_test_nest_lock(omp_nest_lock_t*) : является неблокирующей версией API-функции omp_set_nest_lock .
• omp_unset_nest_lock(omp_nest_lock_t*) : освобождает блокировку при достижении внутренним счетчиком блокировок значения 0 . В остальных случаях каждый вызов этой API-функции уменьшает значение счетчика.

Детальный контроль выполнения заданий

Мы уже видели, что блок кода, который следует за директивой pragma omp parallel , обрабатывается параллельно всеми потоками. Код внутри этих блоков можно также разделить на категории, которые будут выполняться в указанных потоках. Рассмотрим код в листинге 12.

Листинг 12. Использование прагмы parallel sections

int main() { #pragma omp parallel { cout << "Это выполняется во всех потоках\n"; #pragma omp sections { #pragma omp section { cout << "Это выполняется параллельно\n"; } #pragma omp section { cout << "Последовательный оператор 1\n"; cout << "Это всегда выполняется после оператора 1\n"; } #pragma omp section { cout << "Это тоже выполняется параллельно\n"; } } } }

Код, предшествующий директиве pragma omp sections и следующий сразу за директивой pragma omp parallel , обрабатывается параллельно всеми потоками. При помощи директивы pragma omp sections код, следующий за ней, разбивается на отдельные подсекции. Каждый блок pragma omp section может выполняться отдельным потоком. Тем не менее отдельные операторы внутри блока section всегда выполняются последовательно. В листинге 13 показаны результаты выполнения кода из листинга 12.

Листинг 13. Результаты выполнения кода из листинга 12

tintin$ ./a.out Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется во всех потоках Это выполняется параллельно Последовательный оператор 1 Это тоже выполняется параллельно Это всегда выполняется после оператора 1

В листинге 13 мы снова видим 8 изначально созданных потоков. Для обработки блока pragma omp sections достаточно трех потоков из восьми. Внутри второй секции мы указали порядок, в котором выполняются операторы вывода текста. В этом и заключается смысл использования прагмы sections . При необходимости можно указывать порядок выполнения блоков кода.

Директивы firstprivate и lastprivate в связке с параллельными циклами

В этой статье я уже показывал, как объявить локальную память потока с помощью директивы private . Но как инициализировать локальные переменные потока? Может быть, синхронизировать их со значением переменной главного потока перед дальнейшими действиями? В таких случаях нам пригодится директива firstprivate .

Директива firstprivate

Используя директиву firstprivate(переменная) , можно инициализировать переменную в потоке, присвоив ей любое значение, которое она имела в главном потоке. Рассмотрим код из листинга 14.

Листинг 14. Использование локальной переменной потока, которая не синхронизирована с главным потоком

#include #include int main() { int idx = 100; #pragma omp parallel private(idx) { printf("В потоке %d idx = %d\n", omp_get_thread_num(), idx); } }

Вот что я получил в результате (ваши результаты могут быть другими).

В потоке 1 idx = 1 В потоке 5 idx = 1 В потоке 6 idx = 1 В потоке 0 idx = 0 В потоке 4 idx = 1 В потоке 7 idx = 1 В потоке 2 idx = 1 В потоке 3 idx = 1

В листинге 15 содержится код с использованием директивы firstprivate . Как и ожидалось, вывод результатов показывает, что переменная idx имеет значение 100 во всех потоках.

Листинг 15. Использование директивы firstprivate для инициализации локальных переменных потоков

#include #include int main() { int idx = 100; #pragma omp parallel firstprivate(idx) { printf("В потоке %d idx = %d\n", omp_get_thread_num(), idx); } }

Также обратите внимание на то, что для доступа к идентификатору потока был использован метод omp_get_thread_num() . Этот идентификатор отличается от идентификатора, выводимого командой top операционной системы Linux®, и эта схема является лишь способом, позволяющим OpenMP отслеживать количество потоков. Если вы планируете использовать директиву firstprivate в коде C++ , то обратите внимание на другую ее особенность: переменная, используемая директивой firstprivate , является копирующим конструктором для инициализации самой себя из переменной главного потока, поэтому если копирующий конструктор является приватным для вашего класса, это может привести к неприятным последствиям. Теперь перейдем к рассмотрению директивы lastprivate , которая во многом является обратной стороной монеты.

Директива lastprivate

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

Листинг 16. Параллельный цикл for без синхронизации данных с главным потоком

#include #include int main() { int idx = 100; int main_var = 2120; #pragma omp parallel for private(idx) for (idx = 0; idx < 12; ++idx) { main_var = idx * idx; printf("В потоке %d idx = %d main_var = %d\n", omp_get_thread_num(), idx, main_var); } printf("Возврат в главный поток со значением переменной main_var = %d\n", main_var); }

На моем компьютере с 8 ядрами OpenMP создает для блока parallel for шесть потоков. Каждый поток, в свою очередь, насчитывает по две итерации в цикле. Итоговое значение переменной main_var зависит от потока, который выполнился последним и, следовательно, от значения переменной idx в этом потоке. Другими словами, значение переменной main_var не зависит от последнего значения переменной idx , но зависит от значения, которое содержала переменная idx того потока, который был выполнен в последнюю очередь. Этот пример продемонстрирован в листинге 17.

Листинг 17. Зависимость значения переменной main_var от последнего выполненного потока

В потоке 4 idx = 8 main_var = 64 В потоке 2 idx = 4 main_var = 16 В потоке 5 idx = 10 main_var = 100 В потоке 3 idx = 6 main_var = 36 В потоке 0 idx = 0 main_var = 0 В потоке 1 idx = 2 main_var = 4 В потоке 4 idx = 9 main_var = 81 В потоке 2 idx = 5 main_var = 25 В потоке 5 idx = 11 main_var = 121 В потоке 3 idx = 7 main_var = 49 В потоке 0 idx = 1 main_var = 1 В потоке 1 idx = 3 main_var = 9 Возврат в главный поток со значением переменной main_var = 9

Запустите код из листинга 17 несколько раз, чтобы убедиться в том, что значение переменной main_var в главном потоке всегда зависит от значения переменной idx в последнем выполненном потоке. А что делать, если необходимо синхронизировать значение переменной главного потока с итоговым значением переменной idx в цикле? Здесь на помощь приходит директива lastprivate , работа которой продемонстрирована в листинге 18. Как и в предыдущем случае, запустите код из листинга 18 несколько раз, и вы увидите, что итоговое значение переменной main_var в главном потоке равно 121 (т.е. значению переменной idx в последней итерации цикла).

Листинг 18. Синхронизация с помощью директивы lastprivate

#include #include int main() { int idx = 100; int main_var = 2120; #pragma omp parallel for private(idx) lastprivate(main_var) for (idx = 0; idx < 12; ++idx) { main_var = idx * idx; printf("В потоке %d idx = %d main_var = %d\n", omp_get_thread_num(), idx, main_var); } printf("Возврат в главный поток со значением переменной main_var = %d\n", main_var); }

В листинге 19 представлены результаты выполнения кода из листинга 18.

Листинг 19. Результаты выполнения кода из листинга 18 (обратите внимание на то, что значение main_var always всегда равно 121 в главном потоке)

В потоке 3 idx = 6 main_var = 36 В потоке 2 idx = 4 main_var = 16 В потоке 1 idx = 2 main_var = 4 В потоке 4 idx = 8 main_var = 64 В потоке 5 idx = 10 main_var = 100 В потоке 3 idx = 7 main_var = 49 В потоке 0 idx = 0 main_var = 0 В потоке 2 idx = 5 main_var = 25 В потоке 1 idx = 3 main_var = 9 В потоке 4 idx = 9 main_var = 81 В потоке 5 idx = 11 main_var = 121 В потоке 0 idx = 1 main_var = 1 Возврат в главный поток со значением переменной main_var = 121

И последнее замечание: для поддержки оператора lastprivate объектом C++ требуется, чтобы соответствующий класс содержал доступный публичный метод operator= .

Сортировка слиянием в OpenMP

Рассмотрим реальный пример, в котором OpenMP сокращает время выполнения задачи. Возьмем не слишком оптимизированную версию процедуры merge sort , достаточную для демонстрации преимуществ от использования OpenMP. Этот пример приведен в листинге 20.

Листинг 20. Сортировка слиянием в OpenMP

#include #include #include using namespace std; vector merge(const vector& left, const vector& right) { vector result; unsigned left_it = 0, right_it = 0; while(left_it < left.size() && right_it < right.size()) { if(left < right) { result.push_back(left); left_it++; } else { result.push_back(right); right_it++; } } // Занесение оставшихся данных из обоих векторов в результирующий while(left_it < left.size()) { result.push_back(left); left_it++; } while(right_it < right.size()) { result.push_back(right); right_it++; } return result; } vector mergesort(vector& vec, int threads) { // Условие завершения: список полностью отсортирован, // если он содержит только один элемент. if(vec.size() == 1) { return vec; } // Определяем местоположение среднего элемента в векторе std::vector::iterator middle = vec.begin() + (vec.size() / 2); vector left(vec.begin(), middle); vector right(middle, vec.end()); // Выполнение сортировки слиянием над двумя меньшими векторами if (threads > 1) { #pragma omp parallel sections { #pragma omp section { left = mergesort(left, threads/2); } #pragma omp section { right = mergesort(right, threads - threads/2); } } } else { left = mergesort(left, 1); right = mergesort(right, 1); } return merge(left, right); } int main() { vector v(1000000); for (long i=0; i<1000000; ++i) v[i] = (i * i) % 1000000; v = mergesort(v, 1); for (long i=0; i<1000000; ++i) cout << v[i] << "\n"; }

При использовании 8 потоков для выполнения процедуры merge sort время выполнения сократилось с 3,7 (при использовании одного потока) до 2,1 секунды. Здесь необходимо лишь соблюдать осторожность с количеством потоков. Я начал с 8 потоков, но реальная отдача от их использования может варьироваться в зависимости от конфигурации системы. Если количество потоков не ограничено явным образом, то могут быть созданы сотни, если не тысячи потоков, что с очень высокой вероятностью приведет к снижению производительности системы. Кроме того, при работе с кодом merge sort полезно использовать прагму sections , которая была рассмотрена ранее.

Заключение

На этом я заканчиваю статью. Мы в достаточной мере рассмотрели прагмы parallel и способы создания потоков, убедились в том, что OpenMP уменьшает время выполнения задач и позволяет выполнять синхронизацию и гибкий контроль, а также рассмотрели практический пример с использованием merge sort . Конечно, вам еще предстоит изучить очень многое, и лучше всего поможет в этом Web-сайт проекта OpenMP. Всю дополнительную информацию вы можете найти в разделе .

OpenMP - стандарт программного интерфейса приложений для параллельных систем с общей памятью. Поддерживает языки C, C++, Фортран .

Модель программы в OpenMP

Модель параллельной программы в OpenMP можно сформулировать следующим образом:

• Программа состоит из последовательных и параллельных секций (рис. 2.1).
• В начальный момент времени создается главная нить, выполняющая последовательные секции программы.
• При входе в параллельную секцию выполняется операция fork ,порождающая семейство нитей. Каждая нить имеет свой уникальный числовой идентификатор (главной нити соответствует 0). При распараллеливании циклов все параллельные нити исполняют один код. В общем случае нити могут исполнять различные фрагменты кода.
• При выходе из параллельной секции выполняется операция join .Завершается выполнение всех нитей, кроме главной.

OpenMP составляют следующие компоненты:

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

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

Привязка к языку Fortran

В программах на языке Fortran директивы компилятора , имена подпрограмм и переменных окружения начинаются с OMP . Формат директивы компилятора :

{!|C|*}$OMP директива [оператор_1[, оператор_2, ...]]

Директива начинается в первой (фиксированный формат записи текста языка Fortran 77 ) или произвольной (свободный формат) позиции строки. Допускается продолжение директивы в следующей строке, в этом случае действует стандартное в данной версии языка правило для обозначения строки продолжения (непробельный символ в шестой позиции для фиксированного формата записи и амперсанд для свободного формата).

Пример программы на языке Fortran с использованием OpenMP

program omp_example integer i, k, N real*4 sum, h, x print *, "Please, type in N:" read *, N h = 1.0 / N sum = 0.0 C$OMP PARALLEL DO SCHEDULE(STATIC) REDUCTION(+:sum) do i = 1, N x = i * h sum = sum + 1.e0 * h / (1.e0 + x**2) enddo print *, 4.0 * sum end