GPU编程之数组处理

以下几点是需要明白的：

一：目前，GPU基于数据的并行的处理办法是将块与线程做一个组合，使得每一个线程处理数组中的一个元素！

二：GPU编程是一个数据由主存到显存，经过处理后再把数据送回主存的一个过程！

本人显卡每维包含线程块数量为（65536,65536,65536），相应的每维包含线程数为（1024,1024,64），故可得知线程格的每一维可开启的线程块最大数均为65536，相应线程的最大数为（1024,1024,64）。

                        首先写一个简单的函数来看看GPU处理程序的大致流程！

 

   程序虽然简单但是却说明了开发的大致流程，不管多复杂的算法大概都是按照这个思路来的。

                                            GPU处理一维数组

1）以下为对一维数组并行运算的例子，使用了3种基本方法

#include <cuda_runtime.h>  
#include <iostream>  

#include "book.h"    //该头文件定义了HANDLE_ERROR函数
#include <stdio.h>
using namespace std;

int mode = 0;

#define N 1024

void select()
{
	cout << "select the calculation mode:\n"
		<< "1 -> only by blockIdx\n"
		<< "2 -> only by threadIdx\n"
		<< "3 -> by blockIdx and threadIdx \n";
	cout << "mode =  ";
	cin >> mode;
	cout << endl;
}

__global__ void add_only_blockIdx(int *a, int *b, int *c)       //__global表示该函数可在主机调用，在器件执行；
{                                                 // 另外__device__表示该函数在器件调用，在期间执行
	int idx = blockIdx.x;     //计算位于线程块索引处的数据
	if (idx < N){
		c[idx] = a[idx] + b[idx];
	}
}

__global__ void add_only_threadIdx(int *a, int *b, int *c)
{
	int idx = threadIdx.x;     //计算位于线程索引处的数据
	if (idx < N){
		c[idx] = a[idx] + b[idx];
	}
}

__global__ void add_blockIdx_threadIdx(int *a, int *b, int *c)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;    //将线程块索引与线程索引转为线性
	if (idx < N){
		c[idx] = a[idx] + b[idx];
	}
}

int main()
{
	while (mode< 1 || mode >3){
		select();
	}

	int a[N], b[N], c[N];
	int *dev_a, *dev_b, *dev_c;

	//gpu上分配内存
	HANDLE_ERROR(cudaMalloc((void**)&dev_a, N*sizeof(int)));
	HANDLE_ERROR(cudaMalloc((void**)&dev_b, N*sizeof(int)));
	HANDLE_ERROR(cudaMalloc((void**)&dev_c, N*sizeof(int)));

	//为数组a,b初始化
	for (int i = 0; i < N; ++i){
		a[i] = i;
		b[i] = i;
	}

	//讲数组a,b数据复制至gpu
	(cudaMemcpy(dev_a, a, N*sizeof(int), cudaMemcpyHostToDevice));
	(cudaMemcpy(dev_b, b, N*sizeof(int), cudaMemcpyHostToDevice));

	if (mode == 1)
		add_only_blockIdx << < N, 1 >> >(dev_a, dev_b, dev_c);  //使用N个线程块，其中每个线程块使用1个线程,这里N小于65536，故可行
	else if (mode == 2)
		add_only_threadIdx << < 1, N >> >(dev_a, dev_b, dev_c);  //使用1个线程块，其中该线程块使用1024个线程,这里N小于1024，故可行
	else
		add_blockIdx_threadIdx << < (N + 127) / 128, 128 >> >(dev_a, dev_b, dev_c);  // 这里使用(N+127)/128个线程块，每个线程块128线程，为了
	// 使所 (N+127)/128 * 128 ＞= N         

	//将数组dev_c复制至cpu
	HANDLE_ERROR(cudaMemcpy(c, dev_c, N*sizeof(int), cudaMemcpyDeviceToHost));

	//显示结果
	for (int i = 0; i < N; ++i){
		printf("%d + %d = %d\n", a[i], b[i], c[i]);
	}

	//释放在gpu分配的内存
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_c);
	system("pause");
	return 0;
}

    可以看到，开了1024个线程！但有个小问题是 返回来的数全为0！！ 回头研究下

2）对任意长度的一维数组并行求和运算

       当一维数组长度大于65536*1024= 67108864 时，使用(1)的方法则不可行，因为线程格的每维最大线程块不大于65536，且每维的线程不大于1024（第3维除外），故可使用如下方法。

       思想就是线程的多次利用

#include <cuda_runtime.h>  
#include <iostream>  

#include "book.h"    //该头文件定义了HANDLE_ERROR函数

using namespace std;

#define N 33*1024  //在gpu运算的数组长度受GPU内存的限制，故这里使用33*1024表示长度较大的数组


__global__ void add(int *a, int *b, int *c)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;    //所计算的索引为 128个线程块中的128个线程索引
	// 由于128*128 < 33*1024,故某些线程需执行多次运算
	while (idx < N){
		c[idx] = a[idx] + b[idx];
		idx += gridDim.x * blockDim.x;   //对索引进行递增，递增步长为gridDim.x * blockDim.x，gridDim.x为使</span><span style="font-size:14px;font-family:Arial, Helvetica, sans-serif;">用的线程块总数，</span><span style="font-size:14px;">
	}                                    //blockDim.x为使用的每个线程块中的线程总数
}

int main()
{
	int a[N], b[N], c[N];
	int *dev_a, *dev_b, *dev_c;

	//gpu上分配内存
	HANDLE_ERROR(cudaMalloc((void**)&dev_a, N*sizeof(int)));
	HANDLE_ERROR(cudaMalloc((void**)&dev_b, N*sizeof(int)));
	HANDLE_ERROR(cudaMalloc((void**)&dev_c, N*sizeof(int)));

	//为数组a,b初始化
	for (int i = 0; i < N; ++i){
		a[i] = i;
		b[i] = i;
	}

	//讲数组a,b数据复制至gpu
	(cudaMemcpy(dev_a, a, N*sizeof(int), cudaMemcpyHostToDevice));
	(cudaMemcpy(dev_b, b, N*sizeof(int), cudaMemcpyHostToDevice));


	add << < 128, 128 >> >(dev_a, dev_b, dev_c);  //使用128个线程块，每个线程块使用128个线程,128*128 < 33*1024 


	//将数组dev_c复制至cpu
	HANDLE_ERROR(cudaMemcpy(c, dev_c, N*sizeof(int), cudaMemcpyDeviceToHost));

	//显示结果
	for (int i = 0; i < N; ++i){
		printf("%d + %d = %d\n", a[i], b[i], c[i]);
	}

	//释放在gpu分配的内存
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_c);

	return 0;
}

                                          二维组数（图像）处理并行运算

       对二维图像数据进行并行运算，这里使用经典的Reduce Color例子，其对图像中的每个像素点进行量化，如常见的RGB24图像有256×256×256中颜色，通过Reduce Color将每个通道的像素减少8倍至256/8=32种，则图像只有32×32×32种颜色。假设量化减少的倍数是N，则代码实现时就是简单的value/N*N，通常我们会再加上N/2以得到相邻的N的倍数的中间值！ 图像处理的例子，重点是学习如何操纵二维数组。

//main3.cu
#include <cuda_runtime.h>  
#include <iostream>  

#include "book.h"    //该头文件定义了HANDLE_ERROR函数

#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

     //索引的计算是学习的重点相当于一层一层摞起来了！
__global__ void quantify(uchar *_src_dev, uchar *_dst_dev , int div)
{
	int x_idx = threadIdx.x + blockIdx.x*blockDim.x;           //计算x坐标索引
	int y_idx = threadIdx.y + blockIdx.y*blockDim.y;          //计算y坐标索引
	int idx = x_idx + y_idx * gridDim.x * blockDim.x;        //将(x,y)坐标转换为线性
	_dst_dev[idx ] =  _src_dev[idx] / div*div + div / 2;

}

int main()
{
	Mat src = imread("lena.bmp" );   //读取图像

	Mat dst (src.size(), CV_8UC3);  
	
	uchar *src_data = src.data;          // src_data指针指向图像src的数据
	if (!src.isContinuous())  return -1;   //判断图像数据字节是否填充

	uchar *src_dev, *dst_dev;

	int length = src.rows * src.cols * src.channels();

	//在gpu上分配内存空间
	HANDLE_ERROR(cudaMalloc((void**)&src_dev, length*sizeof(uchar)));
	HANDLE_ERROR(cudaMalloc((void**)&dst_dev, length*sizeof(uchar)));

	HANDLE_ERROR(cudaMemcpy(src_dev, src_data, length*sizeof(uchar), cudaMemcpyHostToDevice));

	dim3 blocks(src.cols*src.channels()/ 32, src.rows/ 32);  //宽和高均为512,所启动线程块与线程数量差距不可过大
	dim3 threads(32, 32);
	quantify << < blocks, threads >> >(src_dev, dst_dev , 64);

	//将数组dst_dev复制至cpu
	HANDLE_ERROR(cudaMemcpy(dst.data, dst_dev, length*sizeof(uchar), cudaMemcpyDeviceToHost));

	cudaFree(src_dev);
	cudaFree(dst_dev);

	imshow("src", src);
	imshow("Dst", dst);
	waitKey(0);


	return 0;
}