Openmp цикл в цикле

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

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

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

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

В статье рассмотрены:

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

Все примеры статьи написаны на С++, использовался компилятор 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 d(a) reduction (+: sum) num_threads(2) { # pragma omp for for (int i = 0; i < n; ++i) sum += a[i]; } return sum; }

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

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

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

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

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

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

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

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

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

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

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

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

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

float rect_integral(const std::<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; }

Читайте также: Bash выйти из цикла

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

рис. 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::<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::<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 d(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. В следующих статьях, возможно, я опишу какие-нибудь другие части этой полезной библиотеки (не забудьте подписаться на рассылку).

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

Вычисление суммы элементов массива [OMP]
Интегрирование методом прямоугольников [OMP]

Ссылки:

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

Источник

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

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

Так что это был мой первый подход:

for(int i=1; i < itemRows; i++){ int itemsIndex = i-1; int itemWeight = integerItems[itemsIndex].weight; int itemWorth = integerItems[itemsIndex].worth; #pragma omp parallel for if(weightColumns > THRESHOLD) for(int c=1; c < weightColumns; c++){ if(c < itemWeight){ table[i][c] = table[i-1][c]; }else{ int worthOfNotUsingItem = table[i-1][c]; int worthOfUsingItem = itemWorth + table[i-1][c-itemWeight]; table[i][c] = worthOfNotUsingItem < worthOfUsingItem ? worthOfUsingItem : worthOfNotUsingItem; } } }

Код работает хорошо, алгоритм решает проблемы правильно.

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

Второй подход:

#pragma omp parallel if(weightColumns > THRESHOLD) { for(int i=1; i < itemRows; i++){ int itemsIndex = i-1; int itemWeight = integerItems[itemsIndex].weight; int itemWorth = integerItems[itemsIndex].worth; #pragma omp for for(int c=1; c < weightColumns; c++){ if(c < itemWeight){ table[i][c] = table[i-1][c]; }else{ int worthOfNotUsingItem = table[i-1][c]; int worthOfUsingItem = itemWorth + table[i-1][c-itemWeight]; table[i][c] = worthOfNotUsingItem < worthOfUsingItem ? worthOfUsingItem : worthOfNotUsingItem; } } } }

Это имеет нежелательный побочный эффект: все внутри параллельного блока теперь будет выполняться n раз, где n – количество доступных ядер. Я уже пытался работать с прагмами а также critical заставить внешний цикл for выполняться в одном потоке, но тогда я не могу рассчитать внутренний цикл по нескольким потокам, если я не открою новый параллельный блок (но тогда не будет никакого увеличения скорости). Но не берите в голову, потому что хорошая вещь: это не влияет на результат. Проблемы по-прежнему решены правильно.

СЕЙЧАС странная вещь: второй подход быстрее, чем первый!

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

Сначала я думал: «Да, возможно, это из-за управления потоками», но потом я прочитал, что OpenMP объединяет созданные экземпляры потоков, что говорит против моего предположения. Затем я отключил оптимизацию компилятора (флаг компилятора -O0), чтобы проверить, имеет ли он какое-либо отношение к этому. Но это не повлияло на измерение.

Кто-нибудь из вас может пролить больше света на это, пожалуйста?

измеренное время для решения задачи о рюкзаке, содержащей 7500 предметов с максимальной емкостью 45000 (создание матрицы 7500×45000, которая намного превосходит используемую переменную THRESHOLD в коде):

Подход 1: ~ 0,88 с
Подход 2: ~ 0.52 с

Заранее спасибо,

phineliner

РЕДАКТИРОВАТЬ:

Измерение более сложной задачи:

К проблеме добавлено 2500 элементов (от 7500 до 10000) (более сложные проблемы в настоящее время не могут быть решены из-за нехватки памяти).

Подход 1: ~ 1,19 с
Подход 2: ~ 0,71 с

EDIT2:

Я ошибся насчет оптимизации компилятора. Это не влияет на измерение. По крайней мере, я не могу воспроизвести разницу, которую я измерял ранее. Я отредактировал текст вопроса в соответствии с этим.

Решение

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

for(int i=1; i<n; i++) { for(int j=0; j<n; j++) { a[i*n+j] += a[(i-1)*n+j]; } }

Один из способов распараллелить это – поменять местами петли

Способ 1:

#pragma omp parallel for for(int j=0; j<n; j++) { for(int i=1; i<n; i++) { a[i*n+j] += a[(i-1)*n+j]; } }

С этим методом каждый поток выполняет все n-1 итерации i внутреннего цикла, но только n/nthreads итерации j, Это эффективно обрабатывает полосы столбцов параллельно. Тем не менее, этот метод очень плохо кешируется.

Другая возможность – распараллелить только внутренний цикл.

Способ 2:

for(int i=1; i<n; i++) { #pragma omp parallel for for(int j=0; j<n; j++) { a[i*n+j] += a[(i-1)*n+j]; } }

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

Другой способ обрабатывать столбцы параллельно, но каждый ряд последовательно:

Способ 3:

#pragma omp parallel for(int i=1; i<n; i++) { #pragma omp for for(int j=0; j<n; j++) { a[i*n+j] += a[(i-1)*n+j]; } }

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

Лучшее решение – это то, которое обрабатывает полосы столбцов параллельно, как метод 1, но все еще дружественно к кешу. Это может быть достигнуто с помощью nowait пункт.

Способ 4:

#pragma omp parallel for(int i=1; i<n; i++) { #pragma omp for nowait for(int j=0; j<n; j++) { a[i*n+j] += a[(i-1)*n+j]; } }

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

Вот время в секундах для n=3000 на моей четырехядерной системе IVB GCC 4.9.2:

method 1: 3.00 method 2: 0.26 method 3: 0.21 method 4: 0.21

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

Из сроков ясно, насколько дружествен к кешу метод 1. Также ясно, что метод 3 быстрее, чем метод 2, и что nowait имеет небольшой эффект в этом случае.

Поскольку метод 2 и метод 3 обрабатывают столбцы в строке параллельно, а строки последовательно, можно ожидать, что их время будет одинаковым. Так почему они отличаются? Позвольте мне сделать несколько замечаний:

Из-за пула потоков потоки не создаются и не уничтожаются для каждой итерации внешнего цикла метода 2, поэтому мне не ясно, каковы дополнительные издержки. Обратите внимание, что OpenMP ничего не говорит о пуле потоков. Это то, что реализует каждый компилятор.
Единственное другое отличие между методом 3 и методом 2 состоит в том, что в методе 2 обрабатывается только главный поток i тогда как в методе 3 каждый поток обрабатывает частный i, Но это кажется мне слишком тривиальным, чтобы объяснить существенную разницу между методами, потому что неявный барьер в методе 3 в любом случае заставляет их синхронизироваться и обрабатывать i это вопрос приращения и условного теста.
Тот факт, что метод 3 не медленнее, чем метод 4, который обрабатывает целые полосы столбцов параллельно, говорит о том, что дополнительные затраты в методе 2 заключаются в уходе и входе в параллельную область для каждой итерации i

Итак, мой вывод заключается в том, что для объяснения того, почему метод 2 намного медленнее, чем метод 3, необходимо изучить реализацию пула потоков. Для GCC, который использует pthreads, это, вероятно, можно объяснить созданием игрушечной модели пула потоков, но у меня пока нет достаточного опыта с этим.

Другие решения

Я думаю, что простая причина в том, что, так как вы размещаете #pragma omp parallel на внешнем уровне области видимости (вторая версия) издержки на вызов потоков менее затратны.

Другими словами, в первой версии вы вызываете создание потоков в первом цикле. itemRows время, тогда как во второй версии вы вызываете создание только один раз.

И я не знаю почему !

Я попытался воспроизвести простой пример, чтобы проиллюстрировать это, используя 4 потока с включенным HT:

#include <iostream> #include <vector> #include <algorithm> #include <omp.h> int main() { std::vector<double> v(10000); std::generate(v.begin(), v.end(), []() { ic double n{0.0}; return n ++;} ); double start = omp_get_w(); #pragma omp parallel // version 2 for (auto& el : v) { double t = el – 1.0; // #pragma omp parallel // version 1 #pragma omp for for (size_t i = 0; i < v.size(); i ++) { el += v[i]; el-= t; } } double end = omp_get_w(); std::cout << ” wall : ” << end – start << std::endl; // for (const auto& el : v) { std::cout << el << “;”; } }

Комментарий / раскомментируйте в соответствии с версией, которую вы хотите. Если вы компилируете с: -std=c++11 -fopenmp -O2 Вы должны увидеть, что версия 2 быстрее.

Демо на Колиру

Живая версия 1 wall : 0.512144

Живая версия 2 wall : 0.333664

Источник