Библиотека OpenMP. Параллельный цикл

Статья ориентирована на тех, кто не знаком с библиотекой OpenMP, но хотел бы познакомиться.

OpenMP — не просто библиотека параллельного программирования, но и стандарт, официально поддерживаемый для языков Си, C++ и Fortran (а неофициально и для других языков, Free Pascal, например [1]). Работает OpenMP только на архитектурах с общей памятью.

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

Не секрет, что большую часть времени выполнения, программы проводят внутри циклов. Наверное поэтому, в OpenMP есть специальная директива параллельного цикла, ей и посвящена большая часть статьи.

В статье рассмотрены 2 примера распараллеливания:

  • вычисление суммы элементов массива;
  • вычисление интеграла методом прямоугольников.

Все примеры статьи написаны на С++, использовался компилятор gcc (но можно использовать и другие, отличаться будут только ключи, передаваемые компилятору). Для поддержки OpenMP, gcc должен принять ключ -fopenmp.

1 Вычисление суммы элементов массива

Как отмечалось выше, OpenMP позволяет распараллеливать нашу программу постепенно, а сначала можно написать последовательную программу и отладить, что мы и сделали:

int sum_arr(int *a, int n) {
  int sum = 0;
  for (int i = 0; i < n; ++i)
    sum += a[i];
  return sum;
}

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

OpenMP нам, также, позволяет распределять работу вручную, но чаще всего это не нужно, ведь есть более удобные средства.

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

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

Поток может использовать как локальные переменные, так и разделяемые. Разделяемые переменные (shared) являются общими для всех потоков, очевидно, с ними надо очень осторожно работать (если хоть один поток изменяет значение такой переменной — все остальные должны ждать — это можно организовать средствами OMP). Все константы являются разделяемыми — в нашем примере, разделяемыми являются переменные «a» и «n».

Поток может содержать набор локальных переменных (опции private и firstprivate), для которых порождаются копии в каждом потоке. Для переменных, объявленных в списке private начальное значение не определено, для firstprivate — берется из главного потока. Все переменные, объявленные внутри параллельной области являются локальными (переменная «i» в нашем примере).

Опция reduction, также, задает локальную переменную (sum), а также, операцию, которая будет выполнена над локальными переменными при выходе из параллельной области («+»). Начальное значение локальных переменных, в этом случае, определяется типом операции (для аддитивных операций — ноль, для мультипликативных — единица).

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

Директива for имеет множество опций, подробнее про которые можно прочитать в толстых учебниках. Внутри цикла можно задать опции private и firstprivate, но кроме того, ряд новых. Например schedule определяет способ распределения итераций между потоками, а nowait — убирает неявную барьерную синхронизацию, которая по умолчанию стоит в конце цикла.

В конце статьи приложен архив с исходным кодом примера. Можно проверить что параллельная программа работает значительно быстрее.

рис. 1 распараллеливание OMP

рис. 1 распараллеливание OMP

Сумма элементов массива считается очень быстро, поэтому чтобы получить результаты, показанные на рисунке, пришлось запустить нашу функцию 10 раз. Большую часть времени работы программы занимает заполнение массива случайными числами, которое выполняется в одном потоке (из за этого результаты чуть смазаны). Однако, не ленивый читатель может без труда распараллелить заполнение массива :).

2 Вычисление интеграла методом прямоугольников

2.1 Задано количество прямоугольников

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

На С++ описанный алгоритм может быть выражен следующим образом (сразу приведен параллельный вариант, т.к. нет ничего нового) :

float rect_integral(const std::function<float(float)> &fun,
                    const float a, const float b, const int n) {
  float h = fabs((b - a) / n);

  float sum = 0;
#pragma omp parallel reduction (+: sum)
{
# pragma omp for
  for (int i = 0; i < n; ++i)
    sum += fun(a + i * h) * h;
} // pragma omp parallel

  return sum;
}

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

рис. 2 деградация точности при большом количестве прямоугольников

рис. 2 деградация точности при большом количестве прямоугольников

В качестве аргумента программа на рис. 2 принимает количество прямоугольников, интегрируется функция x * x на интервале [-1, 1], точное значение интеграла 2/3. Чем больше прямоугольников на ограниченном интервале — тем меньше каждый прямоугольник, следовательно, прямоугольники «плотнее прилегают к графику» и точность должна расти. Точность действительно повышается, мы видим это при увеличении количества прямоугольников с 10 до 100, однако, при очень большом их количестве точность резко падает. OpenMP тут оказывается и не причем (обратите внимание, при компиляции не использовался ключ -fopenmp).

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

OpenMP тут вроде бы и не причем, но оказывается так, что на точность может влиять количество работающих потоков и порядок вычисления. Читатель может убедиться в этом, например последовательно вычислив сумму ряда 1/(x * x) сначала при изменении x от 1 до 100000000, а затем, в обратном порядке. Результаты вычислений будут отличаться, и вычисление в обратном порядке дает более точный результат (если не понятно почему и очень интересно — сообщите, я напишу статью по этой теме). OpenMP не гарантирует определенный порядок вычислений, поэтому и может появляться неожиданная погрешность, на это многократно указывается в некоторых источниках [2].

2.2 Задана точность интегрирования

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

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

float rect_eps_integral(const std::function<float(float)> &fun,
                    const float a, const float b, const float eps) {
  float h = fabs(b - a);
  float sum = h * fun(a);

  while (true) {
    h /= 2; // дробление шага

    float newSum = 0;
    for (float x = a; x < b; x += h)
      newSum += fun(x) * h; // интегрируем с заданным шагом

    if (fabs(sum - newSum) < eps)
      break; // достигнута нужная точность
    sum = newSum;
  }

  return sum;
}

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

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

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

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

Для использования OpenMP в этом примере достаточно определить заранее количество итераций внутреннего цикла и использовать в нем целочисленный счетчик:

float rect_eps_integral_2(const std::function<float(float)> &fun,
                    const float a, const float b, const float eps) {
  int n = 1; // в начале считаем только один прямоугольник
  float sum = fabs(b - a) * fun(a);
  float newSum;

  while (true) {
    n *= 2; // дробление шага
    float h = fabs(b - a) / n;

    newSum = 0;
#pragma omp parallel shared(h) reduction (+: newSum) num_threads(2)
{
# pragma omp for // параллельное интегрирование
    for (int i = 0; i < n; ++i)
      newSum += fun(i * h + a) * h;
} // pragma omp parallel

    if (fabs(sum - newSum) < eps) // проверка точности
      break;
    sum = newSum;
  }

  return sum;
}

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

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

Исходный код:

Ссылки:

  1. Поддержка OpenMP для языка Free Pascal / http://wiki.freepascal.org/OpenMP_support
  2. Антонов А.С. Технологии параллельного программирования MPI и OpenMP: Учеб. пособие. Предисл.: В.А.Садовничий. — М.: Издательство Московского университета, 2012.-344 с.-(Серия «Суперкомпьютерное образование»). ISBN 978-5-211-06343-3.
  3. OpenMP на msdn / https://msdn.microsoft.com/ru-ru/library/dd335940.aspx
  4. официальный сайт OpenMP / http://www.openmp.org/
  5. OpenMP на parallel / http://parallel.ru/tech/tech_dev/openmp.html

2 thoughts on “Библиотека OpenMP. Параллельный цикл

  1. Борис

    Здравствуйте, вы написали статью на счет изменения результата от порядка вычисления в цикле ( конец 2.1.) ?

    Reply

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

*

code