基于GPU的归并算法实现

很久不写程序，犯了很多的错误！

逻辑条理不清晰，分析不够全面 ，蛋疼了两天才弄出来！并且为了，逻辑上的简洁，目前的版本只能用于一个block！

由于CUDA中，并不能实现全局线程的同步，所以，多个block线程的同步没法做到，只能是在每个block完成工作后，重新开启kernel函数，进一步归并，这样就到了更上一层。

虽然没多少人看，还是说一下！

我习惯在代码中附上自己的分析，思路，已经犯过的错误！相信以后看起来会更容易！

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include<stdlib.h>
#include<iostream>
#include<device_functions.h>
#include<string.h>

#include <stdio.h>
using namespace std;

#define N 256//定义有N个数需要排序
#define size N*sizeof(float)
__global__ void mergeSort_kernel(float *d_a,float *d_b);
__device__ void mergeSort(float *d_a,float *d_b,int tid,int stride);
void showData(float *p);

int main()
{
	cudaError_t err;
	err=cudaSuccess;//定义成功码
	float *h_a,*h_b;

	h_a=(float*)malloc(size);
	srand(2014);
	h_b=(float*)malloc(size);memset(h_b,0,size);
	for(int j=0, i=256;i<512;i++,j++)
	{
		//h_a[i]=rand();
		h_a[j]=i;
	}
	cout<<"随机数初始化完成！"<<endl;
	showData(h_a);
	cout<<endl;
	
	float *d_a,*d_b;
	err=cudaMalloc((void**)&d_a,size);
	if(err!=cudaSuccess)
	{
		cout<<"d_a分配空间失败！"<<endl;
		exit(EXIT_FAILURE);
	}
	err=cudaMalloc((void**)&d_b,size);
	if(err!=cudaSuccess)
	{
		cout<<"d_b分配空间失败！"<<endl;
		exit(EXIT_FAILURE);
	}
	err=cudaMemset(d_b,0,size);
	if(err!=cudaSuccess)
	{
		cout<<"d_b初始化失败！"<<endl;
		exit(EXIT_FAILURE);
	}
	err=cudaMemcpy(d_a,h_a,size,cudaMemcpyHostToDevice);
	if(err!=cudaSuccess)
	{
		cout<<"主机拷贝数据到设备失败！"<<endl;
		exit(EXIT_FAILURE);
	}
	
	dim3 blocks((N+255)/255);
	dim3 threads(256);
	cout<<"进入核函数!"<<endl;
	mergeSort_kernel<<<1,threads>>>(d_a,d_b);

	err=cudaMemcpy(h_b,d_b,size,cudaMemcpyDeviceToHost);
	if(err!=cudaSuccess)
	{
		cout<<"数据拷贝回主机失败！"<<endl;
		exit(EXIT_FAILURE);
	}
	//显示结果，检查是否正确
	cout<<"排序后的输出结果："<<endl;
	showData(h_b);
	free(h_a);
	free(h_b);
	cudaFree(d_a);
	cudaFree(d_b);
	return 0;

}
__global__ void mergeSort_kernel(float *d_a,float *d_b)
{
	//第一步：先判断需要循环的次数
	/***********************************
	*	首先，我这里的思路是按照先将所有元素分解为一个，暂时只考虑实现，不考虑线程资源浪费的问题，
	*所以，每次的合并都会让线程数减半，或者还会多一个出来，因为会考虑到元素个数为奇数的情况，这里
	*计算次数的思路可以理解为：最开始元素的跨度stride为1，因为分解为只有一个元素了，第一次合并后，
	*跨度就变成了2，第二次合并后跨度变成了4，第三次合并则跨度变成了8，当合并到最后只有一个线程来
	*处理的时候，这个时候的跨度就会大于或者等于N了。其实，就可以理解为，最后处理的这个线程其实就是
	*第一个线程，而此时要合并的数据已经是整个数组了，所以，本次循环做完就不必再循环了。一个显然的结果
	*就是，5个元素，和8个数据，所需要的合并次数都是3次。
	************************************/
	int stride=1;//跨度，再每次的循环后，跨度会变成两倍
	int count=0;//用来计数这里的N个数据，共需要多少次循环
	for(;stride<N;count++,stride*=2);//此循环做完，就知道了到底需要循环的次数，也就是合并需要进行的次数
	//下面的for循环，就是真正的合并操作了，与线程紧密相关
	int tid=blockIdx.x*blockDim.x+threadIdx.x;//每个线程的id，这里都只是一维的，所以线程id表述起来很简单，重点在于每次执行任务的线程不相同
	stride=1;//前面的stride已经用过了，这里要重新归为1
	for(int i=0,temp=1;i<1;i++)//总共要进行count次合并操作//先只执行一次，看看结果
	{
		//每次筛选线程的时候，都是除以2的num次方，第一次除2，第二次除4，第三次除8……
		temp*=2;//线程数减半，选择线程！每次循环刚好乘2，下面再除的时候，线程就要减半了
		//stride*=2;//线程减半的时候，每合并一次跨度就要乘2，合并后才乘2，这里的跨度实际上指每个（就可以认为是第一个）线程处理数据的长度
		/**********************************************
		*下面函数完成的功能，就是在每次的循环中，先将id除2，这样只有一半的线程是需要工作的，
		*这些线程要完成的任务，就是将上一步中两个线程分别指向的块，排序合并起来
		*注意，是先选择出线程，再在进行合并,排序的操作在合并的时候进行
		***********************************************/
		if(tid%temp==0)//刚好都是线程ID为偶数的线程，因为线程从零开始，第一次循环是tid%2，第二次就是tid%4，第三次就是tid%8
		{
			mergeSort(d_a,d_b,tid,stride);//合并
			stride*=2;//先合并，才再将跨度翻倍
			__syncthreads();
		}
	}
	//当合并操作的循环做完了之后，就表示数据已经排好序了,并且是在d_a中，d_b中也存放着同样的数据
}
//在这个函数里，做的是合并排序操作！具体说，就是将两块区域的内存内容逐一拿出来比较放回去,区间长度由stride决定
//其实比较的两个区间的起始为：tid、tid+stride，区间长度为stride
__device__ void mergeSort(float *d_a,float *d_b,int tid,int stride)
{
	if((tid+stride)<N)//必须是小于N，数组中元素也从0开始排，所以等于N的元素也是不存在的，最大的元素也就是N-1，以5个元素的数组举例就很明显了
						//也就是最后一个线程，如果没有与之配对的区域的话，是不工作的
						//这里本质上也是在限定工作的线程,最后的一个奇数线程不干活，但是值也得拷贝到d_b中
	{
		int m,n,k;//这两个参量作为在比较的时候，两个区域指针的移动标志,k是在d_b中的移动标志
		m=tid;n=tid+stride;k=m;
 		for(int i=0;i<2*stride;i++)//这个循环是为了stride这样一个区间的比较,重点在于下面的交换
									//本质还是在想d_b中插入数据，怎么可能是stride次，按最大的算，应该是2倍的stride
		{
			if(m<(tid+stride)&&n<(tid+2*stride)&&n<N)//小于N表明可以进行比较并交换//这里算是理解了为什么要buffer了！比较交换后的数据往哪里放！我这里就直接放在d_b中，
								//但是，因为我每次的输入数据都是d_a，我这里不考虑效率问题，采取直接再将d_b拷回d_a的策略，否则，就要对d_a，d_b进行交换了
								//除了防止取值越界情况，还有当其中一个区域值已经取完的时候，也应该直接拷贝另一个区域剩下的数据了
			{
				(d_a[m]<=d_a
)?(d_b[k]=d_a[m],m++,k++):(d_b[k]=d_a
,n++,k++);//从小到大的顺序
			}
			else//大于或者等于N则表明，最后一个线程现在所取得值已经越界了，此时直接将剩余数据直接拷入即可,m,n大于上述值了，也表示其中一个区域的取值要越界了
			{
				if(m>=(tid+stride))//前面的区域先排完了的话
				{
					for(;n<(tid+2*stride)&&n<N;n++,k++)//把后面区域中还剩下的直接复制过来,n大于N的时候，就表示越界了
					{
						d_b[k]=d_a
;
					}
					break;//当有一部分提前结束，则提前跳出
				}
				else if(n>=(tid+2*stride)||n>=N)//后面区域先排完了，把前面区域的直接赋值过去//如果是后面的区域越界了的话，处理方式是一样的
				{
					for(;m<tid+stride;m++,k++)//由于取值范围由m限定住了，所以，不会出现数组大于N的情况
					{
						d_b[k]=d_a[m];
					}
					break;
				}
				
			}
		}
		//当交换做完了之后，我这里为了省事，不考虑效率，再将d_b的值复制到d_a中去，因为下次的循环还是d_a为输入数据，省了事，但是牺牲了效率
		for(int p=0;p<N;p++)
		{
			d_a[p]=d_b[p];
		}
	}
	else//对于前面选择的要干活的线程，比如是偶数号线程，但是tid+stride>=N的是不干活的，直接将值拷贝入d_b
	{
		d_b[tid]=d_a[tid];
	}

}
void showData(float *p)
{
	for(int i=0;i<N;i++)
		cout<<p[i]<<" ";
}