Wielowątkowe Sito Alkina (CUDA)

Question

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <conio.h>
#include <iostream>

#define PrimNums 158
unsigned int primenums[PrimNums] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929 };

#define ThreadsPerBlock 512
#define big long long int
#define uint unsigned int 

const big IntSize  = (sizeof(int) * 8);

using namespace std;

__device__ inline bool pos(int* a, big pos) {
	return (a[pos / IntSize] >> (pos % IntSize)) & 1;
}
__device__ inline void set(int* a, big pos,bool val) {
	if (val) { //ustawia pos'ty bit w ciągu a na wartość val
		atomicOr(&a[pos / IntSize], (1 << (pos % IntSize))); 
	}
	else {
		atomicAnd(&a[pos / IntSize], ~((char)(1 << (pos % IntSize))));
	}
}
inline bool Pos(int* a, big pos) {
	return (a[pos / IntSize] >> (pos % IntSize)) & 1;
}
inline void Set(int* a, big pos, bool val) {
	if (val) {
		a[pos / IntSize] |= (1 << (pos % IntSize));
	}
	else {
		a[pos / IntSize] &= ~((char)(1 << (pos % IntSize)));
	}
}

__global__ void memoryset(int *a, big num, const unsigned int AllocSize)
{
	big tid = blockDim.x*blockIdx.x + threadIdx.x;
	if (tid == 0) {
		a[0] = 38;
	}
	else {
		while (tid < AllocSize) {
			a[tid] = 0;
			tid += blockDim.x*gridDim.x;
		}
	}
}



__global__ void AtkinSieveKernel(int *a, big num, big sqrtnum, const unsigned int AllocSize)
{
	uint x = blockDim.x*blockIdx.x + threadIdx.x+1;
	uint y = blockDim.y*blockIdx.y + threadIdx.y+1;
	const uint mx = blockDim.x*gridDim.x;
	const uint my = blockDim.y*gridDim.y;
	uint mod;
	uint n;
	
	for (uint i = x; i <= sqrtnum; i += mx) {
		for (uint j = y; j <= sqrtnum; j += my) {
			n = 4*i*i + j*j;
			if (n <= num) {
				mod = n % 60;
				if (mod == 1 || mod == 13 || mod == 17 || mod == 29 || mod == 37 || mod == 41 || mod == 49 || mod == 53) {
					set(a, n, true);
				}
			}
			
			n -= i * i;

			if (n <= num) {
				mod = n % 60;
				if (mod == 7 || mod == 19 || mod == 31 || mod == 43) {
					set(a, n, true);
				}
			}
			
			n -= 2*j * j;

			if (n <= num) {
				mod = n % 60;
				if (mod == 11 || mod == 23 || mod == 47 || mod == 59) {
					set(a, n, true);
				}
			}
		}
	}
}

bool SitoAtkina_GPU(const big num) {
	
	const big AllocSize = (num - 1) /IntSize + 1;
	const big sqrtnum = sqrt(num);

	int *a = new int[AllocSize];

	int *gpu_a;
	cudaError Error;
	Error = cudaMalloc((void**)&gpu_a, AllocSize * sizeof(int));
	if (Error != cudaSuccess) {
		cout << "Za malo VRAMU" << endl;
		return false;
	}

	cout << "malloc Size " << AllocSize << endl;

	cudaDeviceProp GPU;
	cudaGetDeviceProperties(&GPU, 0);
	uint blocks = GPU.multiProcessorCount * 2;
	uint threads = ThreadsPerBlock;
	

	memoryset <<< GPU.multiProcessorCount * 2, ThreadsPerBlock >>>(gpu_a, num, AllocSize);
	AtkinSieveKernel <<< blocks, threads >>>(gpu_a, num, sqrtnum, AllocSize);

	cudaMemcpy(a, gpu_a, AllocSize * sizeof(int), cudaMemcpyDeviceToHost);

	for (uint i = 5; i <= sqrtnum; i++)
	{
		if (Pos(a, i))
		{
			for (uint j = i * i; j <= num; j += i)
			{
				Set(a, j, false);
			}
		}
	}
	//output a
	cout << "2 2" << endl << "3 3" << endl;
    uint l = 2;
    for (int i = 5; i < PrimNums; i++)
    {
        if (Pos(a, i)) {
            cout << i << " " << primenums[l] << endl;
            l++;
        }
    }

}

int main(void)
{
	SitoAtkina_GPU(PrimNums);
	getch();
	return 0;
}

Dzień dobry. w ramach praktyki postanowiłem spróbować zaimplementować Sito Alkina w CUDA jeszcze nie go nie optymalizowałem na razie chcę by po prostu działał. i tu pojawia się problem bo moje sito jest dziurawe i na pierwszy rzut oka wszystko jest ok, nawet wygląda jakby działało a tu.. 91 to według sita liczba pierwsza! Nie piszcie jak to jest źle zaimplementowane itp bo jeszcze się tym zajmę (np użycie long long int taaak tak wiem że to nie działa).  
wgl zadaję sobie pytanie czy zrównoleżnienie tego sita jest możliwe/opłacalne? 

Jak już rozmawiamy o sitach to chciałbym się zapytać czy da się temu situ nałożyć dolny limit? tzn np aby wyznaczyło liczby pierwsze od 100mld do 101mld.

Answer 1

Sito z definicji dość trudno zrównoleglić (oczywiście mówię o klasie algorytmów) ze względu na dostęp do wyników działania poprzedniego etapu obliczeń w etapie bieżącym. I jest to także odpowiedź na 2 Twoje pytanie co do granic obliczeń w sicie. Możesz za to łatwo wykonać operacje asynchronicznie i je "rozleniwić". Zanim jednak będziesz zastanawiał nad tą implementacją, proponuję poprawnie zaimplementować samo sito. To że zwraca nieprawidłowy wynik, świadczy o błędzie w algorytmie (na 80%). Napisz testy i dopiero później zrób z tego kernel. Na końcu próbuj optymalizować. Dobre osadzenie obliczeń na GPU to nie jest trywialna sprawa.

PS. Zwróć uwagę na rodzaje typów danych które obsługuje Twoja karta. Możesz się zdziwić także co do przykładów dostarczanych wraz z SDK :)