我想确定有多少 x^2+1 形式的数是素数,因为 1 <= x <= 10^7。我只是想将它与 CUDA 并行化并检查差异,所以我使用了简单的素数检查,我不关心改进它的算法。
我安排了一个网格,并将其滑过我的间隔,将结果记录在每个块的共享内存中,对每个块执行 gpu 缩减,最后执行 cpu 缩减以获得最终结果。
我的问题是,当我更改每个块中的块数和线程数时,输出结果会发生变化。我无法解释的另一件事是,对于 8 个块和每个块 2048 个线程的配置,代码在 100 毫秒内运行,但是当我将线程数减少到 1024 并将块数加倍时,代码将导致超时在从设备到主机的 memcpy 中!!我该如何解释这种行为以及正确性在哪里出现问题?
我正在使用 GTX 480 nvidia gpu。
我的代码是:
#include <stdio.h>
static void HandleError( cudaError_t err, const char *file, int line )
{
if (err != cudaSuccess) {
printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
exit( EXIT_FAILURE );
}
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
#define N 10000000
#define BLOCKS 8
#define THREADS 2048
__device__ int isprime(int x)
{
long long n = (long long)x*x + 1;
for( int p=3; p<=x+1; p+=2 )
if ( n % p == 0 ) return 0;
return 1;
}
__global__ void solve(int n, int* result)
{
__shared__ int ipc[THREADS];
int tid = threadIdx.x;
int x = blockIdx.x*blockDim.x + threadIdx.x + 2;
// sliding grid window over interval of to-be-computed data
int acc = 0;
while( x <= n )
{
if ( isprime(x) ) acc++;
x += blockDim.x*gridDim.x;
}
ipc[tid] = acc;
__syncthreads();
// reduction over each block in parallel
for( int s=blockDim.x/2; s>0; s>>=1 )
{
if ( tid < s )
{
ipc[tid] += ipc[tid+s];
}
__syncthreads();
}
if ( tid == 0 ) result[blockIdx.x] = ipc[0];
}
int main()
{
int *dev;
int res[BLOCKS];
int ans = 0;
HANDLE_ERROR( cudaMalloc((void**)&dev, BLOCKS * sizeof(int)) );
solve<<<BLOCKS, THREADS>>>(N, dev);
HANDLE_ERROR( cudaMemcpy(res, dev, BLOCKS*sizeof(int), cudaMemcpyDeviceToHost) );
// final reduction over results for each block
for( int j=0; j<BLOCKS; j++ )
ans += res[j];
printf("ans = %d\n", ans);
HANDLE_ERROR( cudaFree( dev ) );
return 0;
}