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.

展开阅读全文

没有更多推荐了,返回首页