利用__shfl_up()函数进行数组求和的诡异问题

我定义了一个一维数组,其中有16个int对象,并将它们每4个为一组,在各组中求之前对象的和。例如:
1、2、3、4为前4个对象,所得结果为1、1+2=3、1+2+3=6、1+2+3+4=10,以此类推。
我分别写了两种代码,但结果却不一样,代码1结果错误,代码2结果正确。不同之处在于:代码1中用一个sum变量来存储之前对象的和,然后将sum与本地的对象相加,最后输出;代码2是直接将之前的对象一个个加到本地上,直接改变本地值,最后输出本地值。代码1每组中的第4个结果比代码2的正确结果都要少加一个值(就是每组的第一个对象),按循环的语句来看,确实应该发生这样的事情,我们还是以1、2、3、4为例。1是第0位,所以它不加任何值,还是输出本身1(正确);2是第1位,offset=1时,返回第0位的值1,offset变为12=2时,不返回值,所以输出1+自身2 =3(正确);3是第2位,offset=1时,返回第1位的值2,offset变为12=2时,返回第0位的值1,所以输出2+1+自身3=6(正确);4是第3位,offset=1时,返回第2位的值3,offset变为1*2=2时,返回第1位的值2,所以输出3+2+自身4=9(错误),正确值应该是10,少加第0位的值1;其他组也都是这个问题。那为什么在代码中这些问题就没有了呢?输出的值都是正确的!
我以为这两种方法的结果应该一致,但结果不对。主要是我不理解GPU在执行循环展开和__shfl_up()函数时的运行机制是什么?具体步骤是什么?
请各位大神帮忙看看问题出在哪里了,代码如下:

#include “cuda_runtime.h”
#include “device_launch_parameters.h”

#include <stdio.h>
#include
using namespace std;
template
global void test(int *s_num_pts, int * s_output1, int * s_output2)
{
// The number of warps in a block.
const int NUM_WARPS_PER_BLOCK = NUM_THREADS_PER_BLOCK / warpSize;
// Compute the coordinates of the threads in the block.
const int warp_id = threadIdx.x / warpSize;
const int lane_id = threadIdx.x % warpSize;

//代码1
if( warp_id < 4 )
{
int num_pts = lane_id < NUM_WARPS_PER_BLOCK ? s_num_pts[warp_id4+lane_id] : 0;
int sum = 0;
#pragma unroll
for( int offset = 1 ; offset < NUM_WARPS_PER_BLOCK ; offset
= 2 )
{
int n = __shfl_up( num_pts, offset, NUM_WARPS_PER_BLOCK );
if( lane_id >= offset )
sum += n;
}
if( lane_id < NUM_WARPS_PER_BLOCK )
s_output1[warp_id4+lane_id] = sum + num_pts;
}
//代码2
if( warp_id < 4 )
{
int num_pts = lane_id < NUM_WARPS_PER_BLOCK ? s_num_pts[warp_id
4+lane_id] : 0;
#pragma unroll
for( int offset = 1 ; offset < NUM_WARPS_PER_BLOCK ; offset*= 2 )
{
int n = __shfl_up( num_pts, offset, NUM_WARPS_PER_BLOCK );
if( lane_id >= offset )
num_pts += n;
}
if( lane_id < NUM_WARPS_PER_BLOCK )
s_output2[warp_id4+lane_id] = num_pts;
}
}
int main()
{
cudaSetDevice(0);
cudaDeviceProp properties;
cudaGetDeviceProperties( &properties, 0);
int warp_size = properties.warpSize;
const int NUM_THREADS_PER_BLOCK = 128; // Do not use less than 128 threads.
const int NUM_WARPS_PER_BLOCK = NUM_THREADS_PER_BLOCK / warp_size;
int * data = new int [16];
for(int i = 0; i<16; ++i){
data[i]= i + 1;
}
for(int i = 0; i<16; ++i){
cout<<data[i]<<" ";
}
cout<<endl;
int dataout1 = new int [16];
int dataout2 = new int [16];
int
input;
cudaMalloc(&input, 16
sizeof(int));
cudaMemcpy(input,data,16
sizeof(int),cudaMemcpyHostToDevice);
int* output1;
cudaMalloc(&output1, 16sizeof(int));
int
output2;
cudaMalloc(&output2, 16sizeof(int));
test<NUM_THREADS_PER_BLOCK> <<<1, NUM_THREADS_PER_BLOCK>>>( input,output1,output2 );
cudaMemcpy(dataout1,output1,16
sizeof(int),cudaMemcpyDeviceToHost);
for(int i = 0; i<16; ++i){
cout<<dataout1<<" ";
}
cout<<endl;
cudaMemcpy(dataout2,output2,16*sizeof(int),cudaMemcpyDeviceToHost);
for(int i = 0; i<16; ++i){
cout<<dataout2<<" ";
}
cout<<endl;
[/i][/i]
[i][i][i][i] }

结果如下:[/i][/i][/i][/i]
[attach]3119[/attach]

shuffle相关函数是用来在warp内进行数据快速交换的。

他的第三个参数,是分组的宽度,数据只有在warp的这个分组里才会交换。

针对你的问题,你需要分组宽度是4(4组4个int),
而楼主你写成:__shfl_up(…,…, NUM_WARPS_PER_BLOCK);
虽然你的NUM_WARPS_PER_BLOCK是128个线程除以32线程/warp, 依然是4,但这个只是巧合让你不出错。

请立刻将参数NUM_WARPS_PER_BLOCK改成4,即:
for( int offset = 1 ; offset < 4; offset*= 2 )
{
__shfl_up( num_pts, offset, 4);
}

上文是其一,不小心用opera的CTRL+空格切换输入法的时候发出了,下文将继续说其他方面:

其二,在楼主的写法(1)中。因为实际上上文说过你的循环是:
for( int offset = 1 ; offset < 4; offset*= 2 )
{
int n = __shfl_up( num_pts, offset, 4);
if(…) sum += n;
}
我们用组0里的4个线程为例, (1,2,3,4四个int) 显然sum的值:
第一次循环前:0 0 0 0
第一次循环后: 0 1 2 3 (0+0, 0 + 1, 0 + 2, 0+3)
第二次循环后:0 1 3 5 (0+0, 1 + 0, 2 + 1, 3 + 2)
写入前的加法: 1 3 6 9 (0+1, 1 + 2, 3 + 3, 5 + 4)
也就是你看:
输出[0]号值是:线程0的值
输出[1]号值是:线程0的值+线程1的值
输出[2]号值是: 线程0的值+线程1的值+线程2的值
输出[3]号值是:线程1的值+线程2的值+线程3的值(少了线程0的!)
为何会如此?
因为你的循环是*2的。offset分别是1和2. 这样第四个线程(线程3),只能加上和他分别相差1和相差2的,而相差3的第一个线程的值就无法加到他身上了。

那么如何解决?请考虑使用(2)中的规约例子(2种的每个线程的传输到临近线程的值总是最近一次加法以后的!而不是原始的。所以可以)。
或者,如果非要结果正确的话,请将循环的offset *=2 改成:offset++

这是你的第二个问题。下文继续。

其三,使用lane_id是个好习惯,但错误的使用它则不是。

楼主的线程模型是(128,1,1)的block,
那么自然这里:
const int warp_id = threadIdx.x / warpSize;
const int lane_id = threadIdx.x % warpSize;
正常写法是无问题的。

但是楼主忘记了shuffle里面会再次进行虚拟lane_id指派。你分成了4个组,每个组里面的lane_id被当作从0-3重新看待。所以你原来的0-15的lane_id这里不适用,实际上称为了4个单独的0-3的lane_id。

建议的方案:
将lane_id改成threadIdx.x & 3

这样可以取得和虚拟分组里的id一样的效果。

大致这些问题吧。也许还有更多。不过暂时我只看出这么多。

[

谢谢您的悉心回复,收获很大!
但我对您提的第二点问题还有些疑惑:
您所说的代码2中的规约方法为什么不会出现以下这种多次相加的问题呢?
打个比方,还是以1、2、3、4为例,
1、考虑线程[0]的值,原值为1,它不参与循环,所以值不变为1;
2、考虑线程[1]的值,原值为2,offset=1时,返回n=线程[0]的值=1,所以线程[1]的值变为2+1=3;
3、考虑线程[2]的值,原值为3,这里就会出问题:当offset=1时,返回的n=线程[1]的值,按照您的观点,线程[1]经过最近一次的加法,值已经改变为3,所以返回n=3,线程[2]的值变为3+3=6;到此为止值是正确的,但线程[2]中的循环并没有结束,offset=2时,依然满足条件,即返回n=线程[0]的值,所以线程[2]的值变为3+3+1=7;(错误)
4、考虑线程[3]的值,原值为4,也会出现上述问题,线程[3]+线程[2]+线程[1],但此时线程[2]和线程[1]的值都已经改变了,所以线程[3]的值会变成4+7+3+1=15;这显然是错误的。

但其实程序结果是正确的,我想知道程序到底是如何进行规约的呢?麻烦您啦,谢谢!

楼主您好,shuffle的时候,整个warp的线程们是并行执行的(同时), 提供交换的数据和得到交换的结果是一步完成的,没有中间状态。

所以:
初始状态:1 2 3 4
第一次交换后: 1 3(+1) 5(+2) 7(+3)
第二次交换后:1 3 6(+1) 10(+3)

例如第一次交换时候的线程1,它从线程0得到1, 同时将自己的值2给线程2。
而不是先得到1, 进行加法得到2+1 = 3, 再将3给线程2的。

warp中是同时进行此操作的,一步完成。

噢。。。。。。原来是这样,太感谢您了!
祝您一切顺利!

您客气了。服务您是我们的荣幸。

[

您好,我还有一个问题想咨询您一下:
我现在的GPU是GTX660,计算能力是3.0,我看自2.1以后,cuda就支持函数的递归了,那形如以下形式的递归,我的GPU是否支持呢?
globle function(int *i)
{
int * j;

function<<<1,128>>>(j);

}

之所以有这个疑问是因为现在计算能力为3.5的GPU可以支持动态并行,例如cuda5.0中给得基于3.5的一个四叉树构建的例子中,其结构大致为:
globle function(int *i)
{
int * j;

function<<<4,128>>>(j);

}
void main()
{
int * init;

function<<<1,128>>>(init);

}

意思就是可以自身动态的在GPU中Launch 4个kernel,而不用返回数据,让CPU重新Launch kernel。这对于构建树一类的递归算法确实非常方便,因为不需要把很多的数据在CPU与GPU中来回传递(这个过程非常容易出错)。

请问如果无法使用动态并行技术,那么如何高效的实现递归算法?cuda所支持的递归到底是那种形式的?

LZ您好,您给出的代码段:
globle function(int *i)
{
int * j;

function<<<1,128>>>(j);

}


是不能在您的GPU上执行的。

因为在kernel里面启动kernel是计算能力为3.5的GPU特有的功能,目前计算能力3.5的GPU有,telsa K20/K20X 和GEFORCE TITAN。

您的卡的计算能力是3.0的,不支持此特性。

以及您提到的“自2.1以后,cuda就支持函数的递归了”,这个指的是在kernel内部递归调用__device__函数,因为这个实际上是每个线程自己的行为,__device__函数的规模受到一个线程自身所能获取的计算资源的限制,规模是比较小的。

您的计算能力为3.0的显卡只支持这种__device__函数的递归调用,而无法在device端(GPU端)实现kernel的递归(即kernel启动kernel)。

以及在您现有的硬件上,一般只能是将参数传回host端,然后重新启动kernel,以实现kernel级别的递归调用。

如您所说,这个比较麻烦,而且效率可能也会受到影响,需要您仔细撰写代码,并测试效率情况。

或者您也可以考虑修改算法,如果可能的话。

或者您也可以直接购买计算能力为3.5的GPU,以直接利用这一硬件能力,以节省人力资源。

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

非常感谢您,收获很大!

不客气的,祝您编码顺利~