CUDA: GPU高性能运算

盆友说，1000个元素，还不如CPU来的快，对的，很多情况下，数据量并行度不是特别大的情况下，可能CPU来的更快一些，比较设备与主机之间互相调度操作，是会有额外开销的。有人就问了，一个10000个元素的数组是不是上面提供的idea就解决不了啦？对，一个block人都没怎么多，如何完成！这个情况下有两条路可以选择——

第一，我就用一个组的1000人来干活话，每个人让他干10个元素好了。

这个解决方案，我们需要修改的是内核函数：

[cpp] view plaincopy

__global__ void kernelarr(int *dev_arr)
{
int tid = threadId.x;
if(tid < 1000) // 只用0~999号线程
{
//每个线程处理10个元素，比如0号线程处理0、1001、2001、……9001
for(int i = tid; i<N; i=i+1000)
{
dev_arr[tid] ++;
}
}
}

第二，我多用几个组来干这件事情，比如我用10个组，每个组用1000人。

这个解决方案就稍微复杂了一点，注意只是一点点哦~因为，组内部怎么干活和最原始的做法是一样的，不同之处是，我们调遣了10个组去干这件事情。

首先我们来修改我们的主机函数：

[cpp] view plaincopy

int main()
{
……
kerneladd<<<10, 1000>>>(dev_arr);//我们调遣了10个组，每个组用了1000人
……
}

盆友要问了，10个组每个组1000人，你怎么点兵呢？很简单啊，第1组第3个线程出列，第9组第9个线程出列。每个人用组号和组内的编号定了位置。在线程网络中，blockId.x和threadId.x就是对应的组号和组内编号啦，我必须要这里开始形象点表示这个对应关系，如果这个对应关系是这样子的[blockId.x，threadId.x]，那么我们的数组arr[10000]可以这样分配给这10个组去干活：

(0,0)——arr[0]，(0,1)——arr[1]，……(0,999)——arr[999]

(1,0)——arr[0+1*1000]，(1,1)——arr[1+1*1000]，……(1,999)——arr[999+1*1000]

……

(9,0)——arr[0+9*1000]，(9,1)——arr[1+9*1000]，……(9,999)——arr[999+9*1000]

是不是很有规律呢？对的，用blockId.x和threadId.x可以很好的知道哪个线程干哪个元素，这个元素的下表就是threadId.x + 1000*blockId.x。

这里我想说的是，如果我们哪天糊涂了，画一画这个对应关系的表，也许，就更加清楚的知道我们分配的线程对应的处理那些东西啦。

一维线程网络，就先学这么多了。

二维网络线程

eg2：int arr[32][16]二维的数组自增1。

第一个念头，开个32*16个线程好了哇，万事大吉！好吧。但是，朕现在想用二维线程网络来解决，因为朕觉得一个二维的网络去映射一个二维的数组，朕看的更加明了，看不清楚自己的士兵，如何带兵打仗！

我还是画个映射关系：

一个block中，现在是一个二维的thread网络，如果我用了16*16个线程。

(0,0)，(0,1)，……(0,15)

(1,0)，(1,1)，……(1,15)

……

(15,0)，(15,1)，……(15,15)

呀，现在一个组内的人称呼变了嘛，一维网络中，你走到一个小组里，叫3号出列，就出来一个，你现在只是叫3号，没人会出来！这个场景是这样的，现在你班上有两个人同名的人，你只叫名，他们不知道叫谁，你必须叫完整点，把他们的姓也叫出来。所以，二维网络中的(0,3)就是原来一维网络中的3，二维中的(i,j)就是一维中的(j+i*16)。不管怎么样，一个block里面能处理的线程数量总和还是不变的。

一个grid中，block也可以是二维的，一个block中已经用了16*16的thread了，那我们一共就32*16个元素，我们用2个block就行了。

先给出一个代码清单吧，程序员都喜欢看代码，这段代码是我抄袭的。第一次这么完整的放上代码，因为我觉得这个代码可以让我说明我想说的几个问题：

第一，二维数组和二维指针的联系。

第二，二维线程网络。

第三，cuda的一些内存操作，和返回值的判断。

[cpp] view plaincopy

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#define ROWS 32
#define COLS 16
#define CHECK(res) if(res!=cudaSuccess){exit(-1);}
__global__ void Kerneltest(int **da, unsigned int rows, unsigned int cols)
{
unsigned int row = blockDim.y * blockIdx.y + threadIdx.y;
unsigned int col = blockDim.x * blockIdx.x + threadIdx.x;
if (row < rows && col < cols)
{
da[row][col] = row * cols + col;
}
}
int main(int argc, char **argv)
{
int **da = NULL;
int **ha = NULL;
int *dc = NULL;
int *hc = NULL;
cudaError_t res;
int r, c;
bool is_right = true;
res = cudaMalloc((void**)(&da), ROWS * sizeof(int*)); CHECK(res)
res = cudaMalloc((void**)(&dc), ROWS * COLS * sizeof(int)); CHECK(res)
ha = (int**)malloc(ROWS * sizeof(int*));
hc = (int*)malloc(ROWS * COLS * sizeof(int));
for (r = 0; r < ROWS; r++)
{
ha[r] = dc + r * COLS;
}
res = cudaMemcpy((void*)(da), (void*)(ha), ROWS * sizeof(int*), cudaMemcpyHostToDevice); CHECK(res)
dim3 dimBlock(16, 16);
dim3 dimGrid((COLS + dimBlock.x - 1) / (dimBlock.x), (ROWS + dimBlock.y - 1) / (dimBlock.y));
Kerneltest<<<dimGrid, dimBlock>>>(da, ROWS, COLS);
res = cudaMemcpy((void*)(hc), (void*)(dc), ROWS * COLS * sizeof(int), cudaMemcpyDeviceToHost); CHECK(res)
for (r = 0; r < ROWS; r++)
{
for (c = 0; c < COLS; c++)
{
printf("%4d ", hc[r * COLS + c]);
if (hc[r * COLS + c] != (r * COLS + c))
{
is_right = false;
}
}
printf("\n");
}
printf("the result is %s!\n", is_right ? "right" : "false");
cudaFree((void*)da);
cudaFree((void*)dc);
free(ha);
free(hc);
getchar();
return 0;
}

简要的来学习一下二维网络这个知识点，

dim3 dimBlock(16, 16);

//定义block内的thread二维网络为16*16

dim3 dimGrid((COLS + dimBlock.x - 1) / (dimBlock.x), (ROWS + dimBlock.y - 1) / (dimBlock.y));

//定义grid内的block二维网络为1*2

unsigned int row = blockDim.y * blockIdx.y + threadIdx.y;

//二维数组中的行号

unsigned int col = blockDim.x * blockIdx.x + threadIdx.x;

//二维线程中的列号

三维网络线程

dim3定义了三维的结构，但是，貌似二维之内就能处理很多事情啦，所以，我放弃学习三维。网上看到的不支持三维网络是什么意思呢？先放一放。

给自己充充电

同一块显卡，不管你是二维和三维或一维，其计算能力是固定的。比如一个block能处理1024个线程，那么，一维和二维线程网络是不是处理的线程数一样呢？

回答此问题，先给出网络配置的参数形式——<<<Dg,Db,Ns,S>>>，各个参数含义如下：

Dg：定义整个grid的维度，类型Dim3，但是实际上目前显卡支持两个维度，所以，dim3<<Dg.x, Dg.y, 1>>>第z维度默认只能为1，上面显示出这个最大有65536*65536*1，每行有65536个block，每列有65536个block，整个grid中一共有65536*65536*1个block。

Db：定义了每个block的维度，类型Dim3，比如512*512*64，这个可以定义3维尺寸，但是，这个地方是有讲究了，三个维度的积是有上限的，对于计算能力1.0、1.1的GPU，这个值不能大于768，对于1.2、1.3的不能大于1024，对于我们试一试的这块级别高点的，不能大于1536。这个值可以获取哦——maxThreadsPerBlock

Ns：这个是可选参数，设定最多能动态分配的共享内存大小，比如16k，单不需要是，这个值可以省略或写0。

S：也是可选参数，表示流号，默认为0。流这个概念我们这里不说。

接着，我想解决几个你肯定想问的两个问题，因为我看很多人想我这样的问这个问题：

1 block内的thread我们是都饱和使用吗？

答：不要，一般来说，我们开128或256个线程，二维的话就是16*16。

2 grid内一般用几个block呢？

答：牛人告诉我，一般来说是你的流处理器的4倍以上，这样效率最高。

回答这两个问题的解释，我想抄袭牛人的一段解释，解释的好的东西就要推广呀：

GPU的计算核心是以一定数量的Streaming Processor(SP)组成的处理器阵列，NV称之为Texture Processing Clusters(TPC)，每个TPC中又包含一定数量的Streaming Multi-Processor(SM)，每个SM包含8个SP。SP的主要结构为一个ALU（逻辑运算单元），一个FPU（浮点运算单元）以及一个Register File(寄存器堆)。SM内包含有一个Instruction Unit、一个Constant Memory、一个Texture Memory，8192个Register、一个16KB的Share Memory、8个Stream Processor(SP)和两个Special Function Units（SFU）。（GeForce9300M GS只拥有1个SM） Thread是CUDA模型中最基本的运行单元，执行最基本的程序指令。Block是一组协作Thread，Block内部允许共享存储，每个Block最多包含512个Thread。Grid是一组Block，共享全局内存。Kernel是在GPU上执行的核心程序，每一个Grid对应一个Kernel任务。在程序运行的时候，实际上每32个Thread组成一个Warp，每个 warp 块都包含连续的线程，递增线程 ID 。Warp是MP的基本调度单位，每次运行的时候，由于MP数量不同，所以一个Block内的所有Thread不一定全部同时运行，但是每个Warp内的所有Thread一定同时运行。因此，我们在定义Block Size的时候应使其为Warp Size的整数倍，也就是Block Size应为32的整数倍。理论上Thread越多，就越能弥补单个Thread读取数据的latency ，但是当Thread越多，每个Thread可用的寄存器也就越少，严重的时候甚至能造成Kernel无法启动。因此每个Block最少应包含64个Thread，一般选择128或者256，具体视MP数目而定。一个MP最多可以同时运行768个Thread，但每个MP最多包含8个Block，因此要保持100%利用率，Block数目与其Size有如下几种设定方式： Ø 2 blocks x 384 threads Ø 3 blocks x 256 threads Ø 4 blocks x 192 threads Ø 6 blocks x 128 threads Ø 8 blocks x 96 threads

这些电很重要啊，必须要充！不然，我就很难理解为什么网络线程如何分配的。

6 规约思想和同步概念

扩大点说，并行计算是有一种基本思想的，这个算法能解决很多很常规的问题，而且很实用，比如说累加和累积等——规约思想。对于基础的、重要的，我想有必要系统的学习。

我觉得有必要重新复制下之前写的这篇介绍：

http://www.cnblogs.com/viviman/archive/2012/11/21/2780286.html

并行程序的开发有其不同于单核程序的特殊性，算法是重中之重。根据不同业务设计出不同的并行算法，直接影响到程序的效率。因此，如何设计并行程序的算法，似乎成为并编程的最大难点。观其算法，包括cuda sdk的例子和网上的牛人，给出的一些例子，以矩阵和矢量处理为主，深入点的包括fft和julia等数学公式，再高级一点的算是图形处理方面的例子。学习这些算法的思想，免不了有自己的一点点总结。之前学习过omp编程，结合现在的cuda，我觉得要理解并行编程，首先理解划分和规约这两个概念。也许你的算法学的更加扎实。划分是《算法》里面的一个重要思想，将一个大的问题或任务，分解成小问题小任务，各个击破，最后归并结果；规约是《cuda**》书上介绍的一个入门的重要思想，规约算法(reduction)用来求连加、连乘、最值等，应用广泛。每次循环参加运算的线程减少一半。不管算法的思想如何花样，万变不离其中的一点--将一个大的任务分解成小的任务集合，分解原则是粒度合适尽量小、数据相关性尽量小。如此而已。因为，我们用GPU是为了加速，要加速必须提高执行任务的并行度！明白这个道理，那么我们将绞尽脑汁地去想方设法分析自己手上的任务，分解、分解、分解！这里拿规约来说事情，因为，规约这个东西，似乎可以拿来单做9*9乘法表来熟悉，熟悉了基础的口诀，那么99*99的难题也会迎刃而解。ex：矢量加法，要实现N=64*256长度的矢量的累加和。假设a+b计算一次耗时t。

cpu计算：显然单核的话需要64*256*t。我们容忍不了。

gpu计算：最初的设想，我们如果有个gpu能同时跑N/2个线程，我们这N/2个线程同时跑，那么不是只需要t时间就能将N个数相加编程N/2个数相加了吗？对的。这一轮我们用了t时间；接着的想法，我们不断的递归这个过程，能发现吗？第二轮，我们用N/2/2个线程同时跑，剩下N/2/2个数相加，这一轮我们同样用了t时间；一直这样想下去吧，最后一轮，我们用了1个线程跑，剩下1个数啦，这就是我们的结果！每一轮时间都为t，那么理想情况，我们用了多少轮这样的计算呢？计算次数=log(N)=6*8=48，对的，只用了48轮，也就是说，我们花了48*t的时间！

规约就是这样，很简单，很好用，我们且不管程序后期的优化，单从这个算法分析上来说，从时间复杂度N降到了logN，这在常规算法上，要提高成这样的效率，是不得了的，这是指数级别的效率提高！所以，你会发现，GPU有CPU无法取代的得天独厚的优势——处理单元真心多啊！

规约求和的核函数代码如下：

[cpp] view plaincopy

__global__ void RowSum(float* A, float* B)
{
int bid = blockIdx.x; int tid = threadIdx.x;
__shared__ s_data[128]; //read data to shared memory
s_data[tid] = A[bid*128 + tid];
__synctheads(); //sync
for(int i=64; i>0; i/=2)
{
if(tid<i) s_data[tid] = s_data[tid] + s_data[tid+i] ;
__synctheads();
}
if(tid==0)
B[bid] = s_data[0];
}

这个例子还让我学到另一个东西——同步！我先不说同步是什么，你听我说个故事：我们调遣了10个小组从南京去日本打仗，我们的约定是，10个组可以自己行动，所有组在第三天在上海机场会合，然后一起去日本。这件事情肯定是需要处理的，不能第1组到了上海就先去日本了，这些先到的组，唯一可以做的事情是——等待！这个先来后到的事情，需要统一管理的时候，必须同步一下，在上海这个地方，大家统一下步调，快的组等等慢的组，然后一起干接下去的旅程。

是不是很好理解，这就是同步在生活中的例子，应该这样说，计算机的所有机制和算法很多都是源于生活！结合起来，理解起来会简单一点。

在CUDA中，我们的同步机制用处大吗？又是如何用的呢？我告诉你，一个正常规模的工程中，一般来说数据都会有先来后到的关系，这一个计算结果可能是提供给另一个线程用的，这种依赖关系存在，会造成同步的应用。

__synctheads()这句话的作用是，这个block中的所有线程都执行到此的时候，都听下来，等所有都执行到这个地方的时候，再往下执行。

7 撬开编程的锁

对于存在竞争的数据，CUDA提供了原子操作函数——ATOM操作。

先亮出使用的例子：

[cpp] view plaincopy

__global__ void kernelfun()
{
__shared__ int i=0;
atomicAdd(&i, 1);
}

如果没有加互斥机制，则同一个half warp内的线程将对i的操作混淆林乱。

用原子操作函数，可以很简单的编写自己的锁，SKD中有给出的锁结构体如下：

[cpp] view plaincopy

#ifndef __LOCK_H__
#define __LOCK_H__
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "atomic_functions.h"
struct Lock {
int *mutex;
Lock( void ) {
HANDLE_ERROR( cudaMalloc( (void**)&mutex, sizeof(int) ) );
HANDLE_ERROR( cudaMemset( mutex, 0, sizeof(int) ) );
}
~Lock( void ) {
cudaFree( mutex );
}
__device__ void lock( void ) {
while( atomicCAS( mutex, 0, 1 ) != 0 );
}
__device__ void unlock( void ) {
atomicExch( mutex, 0 );
}
};
#endif

8 CUDA软件体系结构

9 利用好现有的资源

如果连开方运算都需要自己去编写程序实现，那么我相信程序员这个职业将会缩水，没有人愿意去干这种活。我想，程序员需要学会“偷懒”，现有的资源必须学会高效率的使用。当c++出现了STL库，c++程序员的开发效率可以说倍增，而且程序稳定性更高。

CUDA有提供给我们什么了吗？给了，其实给了很多。

先介绍几个库：CUFFT、CUBLAS、CUDPP。

这里我先不详细学习这些库里到底有哪些函数，但是，大方向是需要了解的，不然找都不知道去哪儿找。CUFFT是傅里叶变换的库，CUBLAS提供了基本的矩阵和向量运算，CUDPP提供了常用的并行排序、搜索等。

CUDA4.0以上，提供了一个类似STL的模板库，初步窥探，只是一个类似vector的模板类型。有map吗？map其实是一个散列表，可以用hashtable去实现这项机制。

SDK里面有很多例子，包括一些通用的基本操作，比如InitCUDA等，都可以固化成函数组件，供新程序的调用。

具体的一些可以固化的东西，我将在以后的学习中归纳总结，丰富自己的CUDA库！

http://blog.csdn.net/huangfengxiao/article/details/8732789

http://blog.csdn.net/huangfengxiao/article/details/8732790

http://blog.csdn.net/huangfengxiao/article/details/8732791