Ищем подматрицу с максимальной суммой элементов? - Найдем и распараллелим!

В статье приводится алгоритм решения задачи, решавшейся по условиям конкурса Acceler8, проводившегося с 15 октября по 20 ноября 2011, организатор – компания Intel.

Основная идея конкурса – получить минимальное время работы программы, хорошую масштабируемость (максимальное ускорение) на всех архитектурах от двухядерной машины среднего пользователя до 40-ядерной машины с общей памятью из Manycore Testing Lab (MTL) от Intel.

Итак, в чем же задача? - Поиск максимальной подматрицы

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

Приведем пример – дана прямоугольная матрица:

-8 -48
-7 -10
3 -55
1 64

Очевидно, прямоугольная подматрица с максимальной суммой элементов в ней есть последний столбец, т.е. матрица с координатами [0][2]…[3][2] (первое число – номер начальной строки, нумерация начинается с нуля, второе – номер начального столбца, третье – конечной строки, четвертое – конечного столбца). Собственно этот ответ и требуется от алгоритма в качестве ответа на исходную матрицу.

Что имеем? Описание входных данных

Запускать будем с 2 параметрами:

Первый параметр командной строки - имя входного файла, второй параметр - имя выходного файла.

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

Первые два числа (M и N) - количество строк и столбцов в матрице, соответственно.

Четыре следующих числа - параметры некоторого псевдослучайного алгоритма, с помощью которого генерируется матрица: seed0, a, b и m соответственно,  – тем не менее алгоритм поиска максимума не зависит от способа генерации входной матрицы. Т.е. Мыне будем привязывать алгоритм к каким-либо предположениям о структуре матрицы, связанным с методом генерации, кроме условия, что в ней есть как отрицательные, так и положительные значения элементов (и нули, конечно :)).

Ограничения для входных параметров (по условиям конкурса):

1 ≤ C ≤ 1000

1 ≤ M, N, a, b, m, seed0 ≤ 20000

Что хотим получить? Описание выходного файла

Для каждого теста в выходной файл выводится строка "Case #x: ", за которой следуют 6 целых чисел (через пробел) и символ переноса строки.

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

Последние 2 числа - сумма элементов в найденной подматрице и площадь этой подматрицы соответственно.

Пример вызова программы

./msp input.txt output.txt

Пример входного файла (input.txt):

2

4 3 10 3 4 17

14 31 11 4 5 18

Пример выходного файла (output.txt):

Case #1: 0 2 3 2 17 4

Case #2: 0 4 5 9 18 36

 

Алгоритмы решения задачи

Последовательный алгоритм решения одномерной задачи

1D:

Очевидно, что «грубая сила» не лучший вариант решения задачи, дающий верхнюю оценку принципиального решения задачи.

За основу был взят алгоритм Кадане 1D, решающий задачу для одномерной задачи за O(N), где N – длина одномерного массива:

Алгоритм Кадане 1D

M=0, t=0
i=1
for j=1 to n do
 	t=t+a[j]
 	if t > M then M=t, (x1,x2)=(i,j)
 	if t ≤ 0 then t=0, i=j+1 //обновить накопление
end for
output M, (x1,x2)

Его работа довольно очевидна: в t накапливается частичная сумма и если t становится больше M, то надо обновить M и положение подмассива. Если t становится отрицательными мы обнуляем накопленное значение, т.к. подмассив с максимальной суммой элементов не может начинаться с меньшего подмассива с отрицательной суммой элементов (в случае рандомного заполнения как положительными, так и отрицательными, очевидно, всегда можно выбрать любое положительное число, а не подмассив с отрицательной суммой элементов;в случае только отрицательных чисел это не так). В приведенном алгоритме M=0 в начале, т.е. если сумма подмассива с максимальной суммой отрицательна, то возвращается пустой подмассив и M=0.

2D:

Обобщение алгоритма на 2D не позволяет в прямом смысле его обобщить, сохранив замечательную скорость поиска максимума – мы можем лишь обобщить его использовав быстрый поиск по 2 измерению и «грубую силу» по 1 измерению:

2D версия Алгоритма Кадане

M = 0, (r1,c1) = (0,0), (r2,c2) = (0,0)
for g = 1 to m do
 	for j = 1 to n do p[j] = 0
 	for i = g to m do
 		t = 0, h = 1
 		for j = 1 to n do
 			p[j] = p[j] + a[i] [j]
 			t = t+p[j]
 			if t > M then M = t, (r1,c1) = (g,h) , (r2,c2) = (i, j)
 			if t ≤ 0 then t=0, i = j+1 //обновить накопление
 		end for
 	end for

end for
output M, (r1,c1), (r2,c2)

Здесь по Y измерению (по номерам строк) мы применяем «грубую силу», а по X измерению (по номерам столбцов) применяется одномерный алгоритм Кадане. Т.к. по Y нам необходимо не меньше, чем O(m2) времени, по X алгоритм Кадане имеет скорость O(n), следовательно общая скорость алгоритма O(m2n).

В программной реализации мы применили алгоритм Кадане по строкам, а не по столбцам, как описано в алгоритме, т.е. далее мы под Х имеем в виду количество столбцов, а под Y – количество строк

Параллелизация и улучшения алгоритма:

Распараллеливание представленного обобщения на 2D алгоритма Кадане:

  • Во-первых, т.к. скорость алгоритма квадратична по m нам необходимо ее уменьшить, если это возможно, т.е. если m>n, то надо искать решение для транспонированной матрицы, что для алгоритма означает, что мы меняем местами m и n, а к 2D матрице, хранящейся по строкам в виде одномерного массива, начинаем обращаться «по столбцам»
  • очевидно, что по X измерению все вычисления независимы везде, кроме сравнивания с текущим максимумом и записи нового значения максимума и его положения. Используя технологию OpenMP, мы можем распараллелить вычисления по внешнему циклу
  • Для начальных Х количество итераций внутреннего цикла от применения «грубой силы», максимально, для конечных Х – минимально. Поэтому при распараллеливании в начале на каждый процессор будет приходиться слишком много работы, а в конце они будут простаивать из-за накладных расходов. Поэтому необходимо, чтобы в начале chunk был равен 1, а в конце – больше 1, чтобы сбалансировать нагрузку. Для OpenMP это решается параметром schedule(guided), при котором размер порции уменьшается с некоторого начального значения до величины chunk (по умолчанию chunk=1) пропорционально количеству ещё не распределённых итераций, делённому на количество нитей, выполняющих цикл. (У нас ситуация с загрузкой прямо противоположная – в начале сильная загрузка, в конце слабая, поэтому мы используем этот параметр вместе с обращением цикла (т.е. отрицательным шагом), т.е. начинаем искать с конца по X).
  • При параллельном выполнении возможны ошибки «гонки данных». Введение критической секции (т.е. области кода, куда может входить одновременно не более одной нити – самым простым решением этих ошибок) – при сравнении с текущим максимумом, обновлением его значения и значений координат начала и конца подмассива – значительно замедлит параллельное выполнение, т.к. большую часть времени нити будут ждать друг друга. Чтобы не вводить их мы создаем массив максимумов (и их координат) – «свой» элемент для каждой нити. Это решает проблему гонки данных и снимает проблему синхронизации. Полная синхронизация нитей (и найденных ими «локальных» максимумов) проводится при выходе из основного вычисления и тривиального поиска максимума в массиве максимум каждой нити.
	int xidx, yidx, endxidx, endyidx;
	int Y = M[p]; //Число строк
	int X = N[p]; //Число столбцов
	
	///* Kadans 2D parallel alghoritm
	int N,M;
    int t = 0;
	int s = 0, k, l,j;
    int z,i,x;
	int numthread = omp_get_num_procs();
	int * S, * x1, * x2, * y1, * y2;
	S = (int*) calloc(numthread, sizeof(int)); //(int*) malloc(sizeof(int)*numthread);
	x1 = (int*) calloc(numthread, sizeof(int));
	x2 = (int*) calloc(numthread, sizeof(int));
	y1 = (int*) calloc(numthread, sizeof(int));
	y2 = (int*) calloc(numthread, sizeof(int));
	omp_set_num_threads(numthread);
//т.к. скорость алгоритма О(X^2*Y), надо в качестве Х надо взять min(X,Y), во втором случае надо обращаться к "транспонированной" матрице
	if (X<=Y){
		N = X;
		M = Y;
		#pragma omp parallel shared(N,M,arr,S,x1,x2,y1,y2) private(z,x,i,t,s,j,k,l)  default(none) 
		{
			#pragma omp for schedule(guided) nowait //shared(N,M,arr,S,x1,x2,y1,y2) private(z,x,i,t,s,j,k,l)  default(none) //
			for(z = N-1; z >-1; z--) //обратное уменьшение шага балансирует нагрузку с chunk > 1, где операций меньше, До chunk=1 в начале, где вычислений больше
			{
				int threadIdx = omp_get_thread_num();
				int * sa;
				sa = (int*) calloc(M,sizeof(int));
				for(x = z; x < N; x++)
				{
					t = 0;
					s = 1<<31;
					j = 0;
					k = 0; l = 0;
					for(i = 0; i < M; i++)
					{
						sa[i] = sa[i] + arr[i*N+x]; 
						t = t + sa[i];
						if( t > s)
						{
							s = t;
							k = i;
							l = j;
						}
						if( t < 0 )
						{
							t = 0;
							j = i + 1;
						}
					}
					if( s > S[threadIdx])
					{
						S[threadIdx] = s;
						x1[threadIdx] = z;
						y1[threadIdx] = k;
						x2[threadIdx] = x;
						y2[threadIdx] = l;
					}
					
				}
				free(sa);
			}
		}
		int max_idx = 0;
		int temp = 1<<31;
		for (i=0;i < numthread;i++)
		{if (temp < S[i]) {temp = S[i]; max_idx=i;}}
		if (x2[max_idx] < x1[max_idx])
		{temp = x1[max_idx]; x1[max_idx] = x2[max_idx]; x2[max_idx] = temp;}
		if(y2[max_idx] < y1[max_idx])
		{temp = y1[max_idx]; y1[max_idx] = y2[max_idx]; y2[max_idx] = temp;}

		FILE* file;
		file = fopen(name, "a");
		if (!file) printf("ERROR: file %s doesn't exist", name);
		fprintf(file, "Case #%i, %i, %i, %i, %i, %i, %i n", p+1, y1[max_idx], x1[max_idx], y2[max_idx], x2[max_idx], S[max_idx], (x2[max_idx]-x1[max_idx]+1)*(y2[max_idx]-y1[max_idx]+1));
		fclose(file);
}

Если m>n то кроме указанных выше изменений в переборе элементов одномерного массива, полученного из исходной матрицы задачи, и замены x<->y при выводе результата код ничем идейно не отличается.

Проблемы распараллеливания 2D алгоритма  Кадане:

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

Тайминги:

Использовались предложенные на форуме:

#include <sys/time.h>
int t = curTime();
fprintf(stderr, "Time of execution: %.3lf sec.n",
(double)(curTime() - t) / 1000.0);

static int curTime() {
    struct timeval now_tv;
    gettimeofday( &now_tv, NULL );
    return (int)now_tv.tv_sec * 1000 + (int)now_tv.tv_usec / 1000;
}

Полученные результаты и недостгнутые цели

  • Транспонирование матрицы (если M<N) дает ускорение в 2 раза для наборов содержащих сильно прямоугольные матрицы (M=3*N). 
  • Использование директивы guided дает ускорение в 1.5 раза.
  • Использование массива максимумов для каждой нити дает ускорение в 5 раз. (Среднее число одновременно работающих потоков остается постоянной, равной удвоенному (с включенным HyperThreading) количеству ядер). Для варианта с критической секции и с общей переменной максимума среднее число одновременно работающих потоков равно 0, т.е. нити всегда друг друга ждут.

На тестах для следующего input.txt:

3
4 3 10 3 4 17
14 31 11 4 5 18
1024 1024 10 3 4 17

Intel Core i3-350M, RAM 3Gb: 9.35 секунд, при этом удалось добиться равномерной загрузки всех 4 потоков (с включенным HyperThreading).

Корректность параллельного алгоритма, т.е. правильность выходных данных для средних датасетов (medium.txt) была проверена на Intel MTL и гарантируется (время выполнения ~23 секунд).

Выводы

Преимущества алгоритма Кадане 2D:

  • Элегантность
  • Простота реализации
  • Относительная простота распараллеливания
  • Хорошая масштабируемость по количеству ядер
  • Большое количество литературы

Недостатки алгоритма Кадане 2D:

  • Далеко не лучшая скорость работы O(m2n) по сравнению с минимальной оценкой идеального алгоритма O(mn). (Лучшую чем O(m2n) оценку имеет, например, алгоритм Такаоки, о котором мы только скажем, что для довольно «маленьких» тестов данного конкурса он слишком громоздок)
  • Принципиальная последовательность алгоритма по второму направлению (столбцы в нашей статье).

Литература:

В процессе поиска было найдено много материалов (ссылки ниже) из открытого доступа.

[1] Tadao Takaoka, Efficient Algorithms for the Maximum Subarray Problem by Distance Matrix Multiplication

[2] KALYAN PERUMALLA, PARALLEL ALGORITHMS FOR MAXIMUM SUBSEQUENCE AND MAXIMUM SUBARRAY

[3] Sung Eun Bae, Sequential and Parallel Algorithms for the Generalized Maximum Subarray Problem

[4] Fredrik Bengtsson C. Jingsen Chen, Efficient Algotithms for k Maximum Sums

[5] Mohammad Bashar, AVERAGE CASE ANALYSIS OF ALGORITHMS FOR THE MAXIMUM SUBARRAY PROBLEM

P.S. Успехов в этом увлекательном деле, Let's go parallel!

For more complete information about compiler optimizations, see our Optimization Notice.