OpenCL多线程累加计算

初学OpenCL，写的不好的地方还请见谅。

本次说的是一个累加和的问题，给定一个一维数组A[n]，然后计算B[n]，其中，B[0]=A[0],B[1]=A[1]+B[0],B[2]=A[2]+B[1]...依此类推。

在我们的CPU中计算的时候使用一个循环来计算，程序是这样的：

B[0]=A[0];

for(int i=1;i<n;i++){

B[i]=B[i-1]+A[i];

}

这样的话就是使用CPU来执行串行计算，如果数据比较多的话效率是非常低的，好的方法就是使用OpenCL来执行并行计算，最好想的并行方法是将两个两个累加，然后再将累加好的数据再两个两个的累加，直到最后只剩下一个数据为止，也就是我们想要的结果，累加的次数于CPU累加的次数一样，但是增加了并行度，比CPU的运行效率要高，累加的图示如下：

上面的方法是从树叶到树根的计算方法，计算出来了最后的值，上层的父节点的值依赖于下层的子节点的值。

但是如果我们想要保存每一次计算的前面的累加和，就要使用从树根到树叶边累加边交换的方法，对于每层累加的次数就等于交换的次数。算法图示如下：

对于这样的小于一个workgroupsize*2的计算不需要任何改进，但是若大于一个workgroupsize*2而我们又需要在一个workgroup内计算完成的话，我们的计算方法肯定需要改进，改进方法也是仁者见仁智者见智，我所使用的就是将数组分段的方法。将给的数组分成独立的X段，每段的元素个数是一个workgroup大小，然后一个线程计算多个元素，分别为每段上本线程id对应的元素（事实上上面的算法是一个线程能计算出来两个元素的，所以每段的元素理论上应该是workgroupsize*2，但是会出错，所以将其调至如此）。然后每段进行上述计算[1]，由于上述算法会将结构推后一位（也就是说B[1]=A[0],B[2]=B[1]+A[1]等），所以每段的最后一个结果会被丢弃，所以在计算过程中需要将每一段的最后的结果存储起来，比如存为C数组，那么第二段的所有结果都要加上C存储的第一段的最后一个结果，同理，第三段的所有结果都要加上C存储的第一段结果和第二段结果之和，依此类推[2]。那么我们在进行这种加之前可以先将C先使用[1]算法累加起来，然后就可以比较简单的进行最后一步加了([2]算法)。具体kernel代码如下所示：

#pragma OPENCL_EXTENSION cl_amd_printf:enable

/*

 * [1]算法

 */

__kernel void column_sum(__global float *src,
	__global float *dst,
	__global float *tmpc,
	__global float *temp,
	const int n1,
	const int globalsize)
{
	const int X=get_global_id(0);
	const int n=n1>>2;

	int count=(n+globalsize-1)/globalsize;

	for(int I=0;I<count;I++){
		barrier(CLK_LOCAL_MEM_FENCE);
		temp[2*X]=src[globalsize*I+2*X];
		temp[2*X+1]=src[globalsize*I+2*X+1];
		int offset=1;

		for(int d=(globalsize>n? n:globalsize)>>1;d>0;d>>=1){
			barrier(CLK_LOCAL_MEM_FENCE);
			if(X<d){
				int ai=offset*(2*X+1)-1;
				int bi=offset*(2*X+2)-1;
				temp[bi]+=temp[ai];
			}
			offset*=2;
		}

		if(X==0){
			if(n>=globalsize)tmpc[I]=temp[globalsize-1];
			temp[(globalsize>n? n:globalsize)-1]=0;
		}
		for(int d=1;d<globalsize&&d<n;d*=2){
			offset>>=1;
			barrier(CLK_LOCAL_MEM_FENCE);
			if(X<d){
				int ai=offset*(2*X+1)-1;
				int bi=offset*(2*X+2)-1;
				float t=temp[ai];
				temp[ai]=temp[bi];
				temp[bi]+=t;
			}
		}
		barrier(CLK_LOCAL_MEM_FENCE);
		dst[globalsize*I+2*X]=temp[2*X];
		dst[globalsize*I+2*X+1]=temp[2*X+1];
	}

}

/*

 * [2]算法

 */

__kernel void scan_sum(__global float *sumtmpc,
	__global float *dst,
	const int n1,
	const int globalsize)
{
	const int X=get_global_id(0);
	if(X<(n1>>2))dst[X]+=sumtmpc[X/globalsize];
}

 

使用此多线程计算和传统使用单线程计算的方法效率比对如下：

全部代码请在http://pan.baidu.com/share/link?shareid=417104&uk=2098130029下载。

参考文献：

[1] Efficient Integral Image Computation on the GPU. Berkin Bilgic, Berthold K.P. Horn, Ichiro Masaki.