这是代码:
#include <stdio.h>
__global__ void para(float* f, int len) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
if (i < len ){
// calculate f[i], f in iteration ith
f[i] = i;
}
}
int main(int argc, char ** argv) {
int inputLength=1024;
float * h_f;
float * d_f;
int size = inputLength*sizeof(float);
h_f = (float *) malloc(size);
cudaMalloc((void**)&d_f , size);
cudaMemcpy(d_f, h_f, size, cudaMemcpyHostToDevice);
dim3 DimGrid((inputLength)/256 +1 , 1 , 1);
dim3 DimBlock(256 , 1, 1);
para<<<DimGrid , DimBlock>>>(d_f , inputLength);
cudaThreadSynchronize();
cudaMemcpy(h_f, d_f, size , cudaMemcpyDeviceToHost);
cudaFree(d_f);
// do parallel reduction
int i;
float sum=0;
for(i=0; i<inputLength; i++)
sum+=h_f[i];
printf("%6.4f\n",sum);
free(h_f);
return 0;
}
并行归约部分可以替换为工作 CUDA 并行求和归约(例如这个)。很快我会花时间改变它。
编辑:
以下是使用 CUDA 执行并行归约的代码:
#include <stdio.h>
__global__ void para(float* f, int len) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
if (i < len ){
// calculate f[i], f in iteration ith
f[i] = i;
}
}
__global__ void strideSum(float *f, int len, int strid){
int i = threadIdx.x + blockDim.x * blockIdx.x;
if(i+strid<len){
f[i]=f[i]+f[i+strid];
}
}
#define BLOCKSIZE 256
int main(int argc, char ** argv) {
int inputLength=4096;
float * h_f;
float * d_f;
int size = inputLength*sizeof(float);
h_f = (float *) malloc(size);
cudaMalloc((void**)&d_f , size);
cudaMemcpy(d_f, h_f, size, cudaMemcpyHostToDevice);
dim3 DimGrid((inputLength)/BLOCKSIZE +1 , 1 , 1);
dim3 DimBlock(BLOCKSIZE , 1, 1);
para<<<DimGrid , DimBlock>>>(d_f , inputLength);
cudaThreadSynchronize();
int i;
float sum=0, d_sum=0;
// serial sum on host. YOU CAN SAFELY COMMENT FOLLOWING COPY AND LOOP. intended for sum validity check.
cudaMemcpy(h_f, d_f, size , cudaMemcpyDeviceToHost);
for(i=0; i<inputLength; i++)
sum+=h_f[i];
// parallel reduction on gpu
for(i=inputLength; i>1; i=i/2){
strideSum<<<((i/BLOCKSIZE)+1),BLOCKSIZE>>>(d_f,i,i/2);
cudaThreadSynchronize();
}
cudaMemcpy(&d_sum, d_f, 1*sizeof(float) , cudaMemcpyDeviceToHost);
printf("Host -> %6.4f, Device -> %6.4f\n",sum,d_sum);
cudaFree(d_f);
free(h_f);
return 0;
}