паралелі непаралельними або пошук простих чисел на GPU

Одним чудовим літнім вечором, я в запалі суперечки мав дурість помітити, що можна написати швидко працююче решето Ератосфена на CUDA. N = 1000000000 (дев'ять нулів) як мету. And the legend has begun…
 
Не буду опускатися в подробиці алгоритму, про нього можна почитати, наприклад, тут і відразу покажу код, яким я мав на той момент:
 
 
#include <iostream>
#include <math.h>

using namespace std;

int main()
{
	double number = 1000000000;
	bool* a = new bool[int(number/2)];
	int i,j,result;

	for (i=0; i<number/2; i++)
		a[i] = true;

	for (i=3; i<=floor(sqrt(number)); i+=2)
		if (a[i/2])
			for (j=i*i; j<=number; j+=i*2)
				a[j/2]=false;

	result = 0;
	for (i=0; i<number/2; i++)
		if (a[i]) result++;

	cout << result << endl;

	delete[] a;

	return 0;
}

Однопотоковий трохи оптимізований код, який працює на 14-15 секунд на Core i3 330M і витрачає дуже багато пам'яті. З нього й почнемо.
 
 

Блочне решето

Запорукою успішного розпаралелювання є поділ завдання на незалежні частини, щоб виключити проблеми з пам'яттю, що виникають при паралельному обігу до неї. І тут вельми до речі доводиться так зване блокове решето, яке дозволяє заощадити пам'ять і розділити велике решето на деяке число маленьких. Можливо це стає завдяки одному цікавому властивості решета, а саме: всі складені числа від кореня з N і до N будуть кратні простим з проміжку від 2 до кореня з N. Інакше кажучи, знайшовши прості до кореня з N (назвемо це попередніми решетом) залишиться тільки для кожного блоку знайти елементи з яких починати викреслювати.
 
Перше, що спадає на думку, це перебір дільників і у випадку подільності починати викреслювати з кроком в дільник. Але навіть за приблизними прикидками подібний підхід виявиться дуже неефективним, навіть з урахуванням усіх можливих оптимізацій.
 
Якщо взяти N = 10000, то тільки для знаходження простих (1229 штук) знадобиться трохи менше тридцяти тисяч перевірок, крім того, не варто забувати, що кожен елемент перевіряється як мінімум один раз, тобто це ще десять тисяч перевірок.
 
До слова, даний варіант був реалізований на CUDA і повісив відеокарту GeForce 9600 GT майже на годину.
 
 

Спроба номер два

На цьому все могло й закінчитися. Але рішення прийшло, як завжди несподівано. І виглядало воно наступним чином.
 
Припустимо, що N = 100, число блоків дорівнює десяти, як і розмір частини. У попередньому решеті будуть чотири елементи: 2, 3, 5, 7. Так як розмір попереднього решета дорівнює розміру частини, то починати можна з другого блоку. Починається він з відомого нам числа 11. Прикладом простого числа за яким буде відбуватися викреслення нам послужить трійка. Двійку не враховую через очевидну оптимізації з її участю, але про це пізніше. Попередній блок закінчився на 10, останній викреслений елемент кратний трійці — дев'ять, знаходиться в одному елементі від кінця блоку. Отже, якщо з трійки відняти цей «хвіст», то можна знайти скільки в наступному блоці залишилося до елемента, який треба бидет викреслити — до 12.
 
І так для кожного простого з попереднього решета. Крім того, оптимізації використовувані в «звичайному» решеті застосовні і тут, а саме: економія пам'яті шляхом виключення парних елементів і викреслення, починаючи з числа в квадраті.
 
 

Від теорії, до практики

Залишилося все теоритические напрацювання зібрати разом. Ось що з цього вийшло:
 
 
for(i=0; i<PARTS; i++)
{
	for(j=0; j<PART_SIZE/2; j++)
		part[j] = true;

	for(j=1; j<result1; j++)
	{
		offset = 0;

		if(part1[j]*part1[j]<=(i+1)*PART_SIZE)
		{
			offset = i*PART_SIZE+part1[j]-i*PART_SIZE%part1[j];
			if(offset%2 == 0) offset += part1[j];
		}
		else break;

		for(k=offset; k<=(i+1)*PART_SIZE; k+=part1[j]*2)
			if(part[(k-(i*PART_SIZE)-1)/2]) part[(k-(i*PART_SIZE)-1)/2] = false;
	}

	if(i==0) result = j - 1;

	for(j=0; j<PART_SIZE/2; j++)
		if(part[j]) result++;
}

Додавати сюди код попереднього решета не вважаю за необхідне, бо він ідентичний «звичайному» решету.
 
 Ця версія решета, назвемо її блокової, працює за 8-9 секунд. image
 

Поділ решета

На досягнутому було вирішено не зупинятися, тому швидко була написана версія працююча в два потоки.
 
 Решето в два потоки, результат — 4.7 секунди. image
Але був один фактор, який уповільнював Виконання покрівельних, другий потік, що працює з другою половиною решета, опинявся медленее основного. Змістивши «розріз» решета, вдалося прискорити виконання майже на півсекунди.
 
 Результат після зсуву — 4.3 секунди. image
 

Решето переїжджає

Перевірений на CPU алгоритм, був перенесений на GPU. Дослідним шляхом було встановлено розміри гріду та частини решета, які дозволяли досягти найбільшу продуктивність.
 
 Решето Ератосфена на CUDA, результат — 4.5 секунди. imageimage
Вихідний код доступний на GitHub .
 
На цьому все, дякую всім за увагу!

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.