CUDA中归约运算

看了好多遍,一直感觉其中temp在while循环中结束后是只得到一个值的,怎么还可能有cache[cacheIndex] = temp;这一步呢?这里要理解了下面也就通了。我不知道我理解错了哪个地方,很迫切想知道。可能问题太简单,但我找了网上的解答,没有回答这个问题的。先谢谢了!
#define imin(a,b) (a<b?a:b)
const int N = 1024* 1024;
const int threadsPerBlock = 256;
const int blocksPerGrid = imin( 512, (N+threadsPerBlock-1) / threadsPerBlock );
global void dot( float *a, float *b, float *c ){
shared float cache[threadsPerBlock];//定义共享存储器,这一步很重要
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0;
while (tid < N) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
// 设置cach[cachIndex]的值
cache[cacheIndex] = temp;

// 块同步处理
__syncthreads();
// for reductions, threadsPerBlock must be a power of 2
// because of the following code
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
}
int main( void ) {
float a, b, c, partial_c;
float dev_a, dev_b, dev_partial_c;
// allocate memory on the CPU side
a = (float
)malloc( N
sizeof(float) );
b = (float
)malloc( N
sizeof(float) );
partial_c = (float
)malloc( blocksPerGrid
sizeof(float) );
// allocate the memory on the GPU
cudaMalloc( (void**)&dev_a,Nsizeof(float) );
cudaMalloc( (void**)&dev_b, N
sizeof(float) );
cudaMalloc( (void**)&dev_partial_c, blocksPerGridsizeof(float) ); // fill in the host memory with data
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = i
2;
}
// copy the arrays ‘a’ and ‘b’ to the GPU
cudaMemcpy( dev_a, a, N* sizeof(float), cudaMemcpyHostToDevice );
cudaMemcpy( dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice ); dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b, dev_partial_c ); // copy the array ‘c’ back from the GPU to the CPU
// finish up on the CPU side

cudaMemcpy( partial_c, dev_partial_c,blocksPerGridsizeof(float),cudaMemcpyDeviceToHost );
c = 0;
for (int i=0; i<blocksPerGrid; i++) {
c += partial_c[i];
}
#define sum_squares(x) (x
(x+1)(2x+1)/6)
printf(“Does GPU value %.6g = %.6g?\n”, c, 2 * sum_squares( (float)(N - 1) ) );
// free memory on the GPU side
cudaFree( dev_a );
cudaFree( dev_b );
cudaFree( dev_partial_c );
// free memory on the CPU side
free( a );
free( b );
free( partial_c );
getchar();
}

LZ您好:

对于每个线程而言,如您所说temp在while循环之后是一个固定值的。

但是您需要考虑到对于一个线程block而言,这里面是多个线程的,每个线程有自己的值,所以这里按照线程编号(cacheIndex = threadIdx.x)将每个线程之前循环的最终结果写到shared memory的数组中并同步,然后进行block多内线程互相协作的规约操作,最终一个block将得到一个唯一的结果,并将此结果写入到global memory中对应的位置。

时刻考虑到您的代码是多个线程在执行的,这也是CUDA多线程编程的一个特点。

大致如此,祝您编码顺利~

版主,对于temp,可以理解为对于不同的线程块、不同的线程都有一个属于自己的temp值吗?是这么理解的不?如果是这样理解的话,那为什么不是个数组,如temp[threadidx.x]?

版主,对于temp,可以理解为对于不同的线程块、不同的线程都有一个属于自己的temp值吗?是这么理解的不?如果是这样理解的话,那为什么不是个数组,如temp[threadidx.x]?

LZ您好:

1:temp变量是您在kernel中定义的,而kernel代码是每个线程都跑一份的,所以这里的temp变量是每个线程私有的,也是每个线程有自己的一个temp。

2:至于为何不是数组,因为您定义的就是一个变量而不是数组。以及,如果您的问题是从block或者grid的角度看为何不是数组,那实际上可以从不同的角度理解:

a)前面说过,temp本身是每个线程独有的变量,那么不会出现a线程访问b线程的temp,所以不会有“temp[threadidx.x]”这样的形式。这样的形式必须是多个线程共享的存储空间的数据才可以,比如global memory和shared memory。

b)根据a)中的叙述,temp的值并不需要维护多个线程的可见性,那么如果维护该变量的存储空间和生命周期就有了更大的灵活性。虽然原则上说,在显存中开辟一块存储空间,以数组的形式保存temp,并且将该数组维护为全部线程的生命周期,这样逻辑上也可以完成之前的需求,但是这并非最佳做法。
事实上,temp是存储于寄存器中的,速度最快,而且对于单个线程而言,当temp不再被使用的时候,可以立即复用该寄存器空间。

大致如此,祝您编码顺利~

版主,首先很感谢您的耐心回答,很受益。我好想懂点了,按您所说的话:当N很大的时候,有tid += blockDim.x * gridDim.x;这个语句,目的就是把N组数据处理完。这样的话,网格内的线程要做多于一次的运算,当第二次全网格运算的时候,temp也不会被覆盖?这时候一个线程保存两个temp值,以此类推。。。
(1)我的理解对不?我表述的很笨拙,不知道您能看懂不。(2)如果是这样的话,那么在寄存器中,它自己也会自己标记出来吗?

LZ您好:

您说的“按您所说的话:当N很大的时候,有tid += blockDim.x * gridDim.x;这个语句,目的就是把N组数据处理完。这样的话,网格内的线程要做多于一次的运算,当第二次全网格运算的时候,temp也不会被覆盖?这时候一个线程保存两个temp值,以此类推。。。”

我没有表达过这个意思,请您重新理解5#相关内容。

顺便说一下,您给出的规约的例子前半段是所有线程都参与的,以每个线程自己累加来实现的规约运算;后半段则是每次参与线程数减半的折半规约。前半段将远大于线程数的数据规约到线程数个数,后半段继续规约到一个确定的结果。

大致如此,希望您能理解。

祝您好运~

谢谢版主。我之前的表达有误,抱歉。对于您说的“您给出的规约的例子前半段是所有线程都参与的,以每个线程自己累加来实现的规约运算”,里面的累加是指的每个线程自己的temp累加吗?对于后半段,我的理解是:对于所有线程开始进行块内的归约。是这样的吧?

LZ您好:

1:是这样的,每个线程将计算结果累加到自己的temp中。
“ while (tid < N)
{
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}”
这段代码清楚地表明了这一点。

2:后一半是块内进行的规约操作。需要指出的是 块与块之间是相对独立的,块内规约只有本块内部的线程参与,其他块内的线程参与其他块的规约操作。块与块之间的执行进度是没有保证的,当一个块执行块内规约的时候,另外一个块可能没开始执行,可能执行到任意位置,也可能已经执行完毕。

所以您“对所有线程开始进行块内规约”这个表述是不确切的,希望您能理解这一点。

祝您好运~

谢谢版主,终于弄懂了,纠结好几天。嘿嘿,再次感谢!

不客气的,欢迎您常来论坛~

祝您编码顺利~