Свёртки в Intel® Cilk Plus

Допустим нам зачем-то нужно найти сумму элементов массива. Мы можем разбить массив на две части, просуммировать каждую часть отдельно и сложить результаты. При этом суммировать эти части можно параллельно. Но суммирование части массива это в точности исходная задача, и каждую часть снова можно разбить на две части и просуммировать каждую часть отдельно, а затем сложить результаты и т. д. Такая стратегия вычислений называется «разделяй и властвуй».

Таким способом можно вычислять много других функций от массивов, ниже в первой части статьи будет приведено математическое объяснение этой идеи, а во второй - как с помощью Intel® Cilk Plus эту идею использовать в своих программах.

Моноиды и гомоморфизмы

Моноидом называется тройка (M, ⋅, e), где M - непустое множество, ⋅ - бинарная ассоциативная операция на множестве M, e - нейтральный по отношению к операции ⋅ элемент множества M.

Так как операция ⋅ является ассоциативной, т. е. для любых x, y и z из M верно, что x ⋅ (y ⋅ z) = (x ⋅ y) ⋅ z, то результат вычисления выражения x1 ⋅ x2 ⋅ … ⋅ xn не зависит от расстановки скобок, поэтому скобки вообще можно не писать.

Примеров моноидов очень много, например, множество целых чисел по сложению (Z, +, 0) является моноидом. Ещё одним примером моноида является списочный моноид. А именно, пусть X - произвольное множество, обозначим символом X совокупность всех конечных последовательностей элементов множества X, включая пустую последовательность [], символом ++ обозначим операцию конкатенации списков. Тогда, как не сложно убедиться, тройка (X, ++, []) является моноидом.

Другие примеры можно найти в статье Моноиды и их приложения: моноидальные вычисления в деревьях.

Далее нам понадобится понятие гомоморфизма. Пусть (M, ⋅MeM) и (N, ⋅NeN) - два произвольных моноида. Функция hM → N называется гомоморфизмом, если h(xMy) = h(x) ⋅Nh(y) для любых x и y из M, и h(eM) = eN.

Заметим, что, нестрого говоря, функция является гомоморфизмом, если её можно протащить через знак бинарной операции. Понятно, что тождественное отображение является гомоморфизмом, но из школьного и университетского курсов математики каждому известно несколько более интересных примеров гомоморфизмов. Например, логарифм произведения равен сумме логарифмов, ln(xy) = ln x + ln y, а в наших терминах это значит, что отображение ln является гомоморфизмом из моноида положительных вещественных чисел по умножению (R+, ⋅, 1) в моноид вещественных чисел по сложению (R, +, 0). Далее, определитель произведения матриц равен произведению определителей этих матриц, det(AB) = det A ⋅ det B, т. е. отображение det является гомоморфизмом из моноида квадратных матриц по умножению (Matn(R), ⋅, I) в моноид вещественных чисел по умножению (R, ⋅, 1). Надеюсь идея ясна. Теперь перейдём к списочным гомоморфизмам.

Любой гомоморфизм h из списочного моноида (X, ++, []) в произвольный моноид (M, ⋅, e) называется списочным гомоморфизмом. Списочные гомоморфизмы обладают следующим полезным свойством. Пусть [x1, x2, …, xn] - произвольный список, тогда h([x1, x2, …, xn]) = h([x1] ++ [x2] ++ … ++ [xn]) = h([x1]) ⋅ h([x2]) ⋅ …⋅ h([xn]). Таким образом, любой списочный гомоморфизм однозначно задаётся своими значениями на одноэлементных списках. Итак, пусть fX → M такая функция, что f(x) = h([x]), тогда однозначно определённый списочный гомоморфизм h будем записывать в виде hom(⋅, fe).

Теперь несколько примеров.

  1. Суммирование элементов списка чисел, произведение элементов списка чисел, нахождение минимального и максимального элементов списка чисел можно рассматривать как списочные гомоморфизмы hom(+, id, 0), hom(⋅, id, 1), hom(max, id, −∞) и hom(min, id, +∞) соответственно. Подобные списочные гомоморфизмы, у которых второй параметр является тождественным отображением id, будем коротко записывать reduce(⋅, e).
  2. Функция length, вычисляющая длину списка, является списочным гомоморфизмом hom(+, one, 0), где one(x) = 1.
  3. Функция, которая применяет некоторую функцию f к каждому элементу списка, тоже является списочным гомоморфизмом, а именно hom(++, g, []), где g(x) = [f(x)]. Такие функции будем коротко записывать как map(f).
  4. Функцию сортировки тоже можно представить в виде списочного гомоморфизма. А именно, пусть merge - функция слияния двух упорядоченных списков, list(x) = [x] - функция, превращающая элемент в одноэлементный список. Тогда гомоморфизм hom(merge, list, []) будет искомым списочным гомоморфизмом.

 

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

Первая теорема о гомоморфизме. Любой списочный гомоморфизм можно представить в виде композиции hom(⋅, fe) = reduce(⋅, e) ∘ map(f).

Именно на этом простом наблюдении основана идея MapReduce. Несложное введение в MapReduce можно найти в статьях MapReduce или подсчеты за пределами возможностей памяти и процессора (попробую без зауми) и MapReduce: более продвинутые примеры, попробуем без зауми.

 

Пусть [x1, x2, …, xn] - произвольный список, левой свёрткой называется функция foldl(⋅LeL)([x1x2, …, xn]) = (…((eLLx1) ⋅Lx2) ⋅L … ⋅Lxn), а правой свёрткой - функция foldr(⋅ReR)([x1x2, …, xn]) =(x1R (x2R … ⋅R (xnReR)…)). Вторая теорема о гомоморфизме утверждает следующее.

 

Вторая теорема о гомоморфизме. Любой списочный гомоморфизм выражается в виде левой и в виде правой свёрток, т. е. hom(⋅, fe) = foldl(⋅Le) = foldr(⋅Re), где xLyx ⋅ f(y), xRyf(x) ⋅ y.

Эта теорема указывает два последовательных способа вычисления списочных гомоморфизмов - в виде левой и в виде правой свёрток. Как называются функции, вычисляющие эти свёртки, в конкретных языках программирования можно найти в статье википедии Fold (higher-order function). Более подробную информацию про свёртки можно найти в статье Евгения Кирпичёва Элементы функциональных языков.

Для нас более интересным является способ параллельного вычисления списочного гомоморфизма h = hom(⋅, f, e). Вспомним, что по определению списочного гомоморфизма h(x ++ y) = h(x) ⋅ h(y), другими словами, можно независимо, а, следовательно, и параллельно, вычислить h(x) и h(y), а затем уже вычислить окончательное значение h(x ++ y). Но каждый из списков x и y мы можем также разбить на две части, и, применяя определение гомоморфизма, вычислить значения h для них независимо, а потом собрать результат, и т. д. Это не что иное, как стратегия «разделяй и властвуй».

Предположим, что операция ⋅ вычисляется за константное время C, а список x состоит из n элементов. Для вычисления левой свёртки нужно в точности n раз вычислить операцию ⋅, т. е. потребуется Cn времени. Если же у нас есть p процессоров, то разобъём список x на p равных частей xi и с помощью левой свёртки для каждого xi вычислим h за время Cn/p, а затем за время Clog2p вычислим окончательный результат. Таким образом, суммарное время равно C(n/p + log2p).

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

Третья теорема о гомоморфизме. Если некоторая функция из списочного моноида в некоторый моноид выражается в виде левой foldl(⋅Le) и в виде правой свёрток foldr(⋅Re), то она является списочным гомоморфизмом.
Доказательства этих трёх теорем можно найти в статье Джереми Гиббонса The Third Homomorphism Theorem, а мы переходим к Intel Cilk Plus.

Свёртки в Intel Cilk Plus

Intel Cilk Plus - это расширение языка C++ для написания многопоточных программ, которое помимо прочего включает в себя следующее:

  1. ключевые слова cilk_spawn и cilk_sync (первое указывает, что функция может выполняться параллельно, второе используется для синхронизации),
  2. ключевое слово cilk_for для распараллеливания циклов,
  3. гиперобъекты для потокобезопасного разделения данных,
  4. специальный синтаксис для записи массивов и операций над массивами, например, поэлементное сложение массивов запишется в виде c[:] = a[:] + b[:].

Здесь будут затронуты только второй и третий пункты. Документацию и другие материалы об Intel Cilk Plus можно найти на сайте /en-us/articles/intel-cilk-plus/. Итак, предположим, что нам нужно подсчитать сумму элементов массива целых чисел. Возможная последовательная реализация приведена ниже.

int sum = 0;
for (int i = 0; i < n; ++i) {
	sum += a[i];
}
std::cout << sum << std::endl;

Заметим, что это есть не что иное, как реализация левой свёртки. Но сумма элементов массива является списочным гомоморфизмом, и, следовательно, допускает также параллельную реализацию, как объяснялось выше. Intel Cilk Plus позволяет сделать такую реализацию безо всяких усилий, достаточно сделать следующее:

  1. подключить заголовочные файлы cilk/reducer_opadd.h и cilk/cilk.h,
  2. изменить тип sum на reducer_opadd<int>,
  3. изменить for на cilk_for,
  4. воспользоваться get_value() для получения результата.
#include <cilk/cilk.h>
#include <cilk/reducer_opadd.h>

...

cilk::reducer_opadd<int> sum;
cilk_for (int i = 0; i < n; ++i) {
	sum += a[i];
}
std::cout << sum.get_value() << std::endl;

При этом произойдёт примерно следующее. Компилятор преобразует тело цикла cilk_forв функцию, реализующую стратегию "разделяй и властвуй", на псевдокоде это запишется следующим образом.

void run_loop(first, last) {
	if ((last - first) < grainsize) {
		for (int i = first; i < last; ++i) {
			LOOP_BODY;
		}
	} else {
		int mid = (last - first) / 2;
		cilk_spawn run_loop(first, mid);
		run_loop(mid, last);
	}
}

А так как переменная sum является редьюсером, то в случае разделения нити на две и их параллельного выполнения, дочерняя будет использовать старое представление, а продолжение получит новое представление, инициализированное нулём, при синхронизации оба преставления будут сложены.
Кроме reducer_opadd Intel Cilk Plus предоставляет большой набор редьюсеров: reducer_list_append, reducer_list_prepend, reducer_min, reducer_max, reducer_min_index, reducer_max_index, reducer_opor, reducer_opand, reducer_opxor, reducer_basic_string, reducer_string, reducer_wstring и reducer_ostream. Если ни один из них не подойдёт, то можно создать свой собственный редьюсер.
Для создания редьюсера нужно создать класс Monoid, реализующий моноид. Этот класс должен содержать следующие функции:

  • reduce(T *left, T *right) вычисляет *left = *left OP *right,
  • identity(T *p) создаёт нейтральный элемент в области памяти, на которую указывает p,
  • destroy(T *p) вызывает деструктор для объекта, на который указывает p,
  • allocate(size) возвращает указатель на область памяти размером size байт,
  • deallocate(p) освобождает память, на которую указывает p.

В большинстве случаев все эти методы реализовывать не нужно, достаточно унаследоваться от класса cilk::monoid_base<T> и реализовать reduce() и, возможно, identity(). Далее, класс, реализующий редьюсер, должен в качестве скрытого члена содержать гиперобъект cilk::reducer<Monoid> imp_и методы для доступа к этому гиперобъекту и для его изменения. Напишем редьюсер для нахождения произведения элементов массива чисел.

#include <cilk/reducer.h>
#include <cilk/cilk.h>
#include <iostream>

class reducer_prod {
public:
	struct Monoid: cilk::monoid_base<double> {
		void reduce(double *left, double *right) const {
			*left *= *right;
		}

		void identity(double* p) const {
			new ((void*) p) double(1.0); 
		}
	};

	reducer_prod(): imp_(1.0) {
	}

	reducer_prod &operator*=(const double &v) {
		imp_.view() *= v;
		return *this;			
	}

    	double get_value() const {
    		return imp_.view();
	}

private:
	cilk::reducer<Monoid> imp_;
};
reducer_prod prod;

cilk_for (int i = 0; i < n; ++i) {
	prod *= a[i];
}

std::cout << prod.get_value() << std::endl;


Для получения подробной информации о возможностях оптимизации компилятора обратитесь к нашему Уведомлению об оптимизации.