《GPU高性能编程CUDA实战 CUDA By Example》

发表于   |   更新于 | 分类于 | 评论数： | 阅读次数：
本文字数： 13k | 阅读时长 ≈ 12 分钟

第一章 为什么需要CUDA

1.1 本章目标

了解历史和发展历程

1.2 并行处理的历史

中央处理器性能的提升逐渐变得困难

1.3 GPU计算的崛起

1.4 CUDA

1.4.1 CUDA架构是什么

cuda架构包含了一个统一的着色器流水线，使得执行通用计算的程序能够对芯片上的每个数学逻辑单元（Arithmetic Logic Unit，ALU）进行排列。另外，NVIDIA在实现ALU时都确保它们满足IEEE单精度浮点数学运算的需求，并且可以使用一个裁剪后的指令集来执行通用计算，而不是仅限于执行图形计算。此外，GPU上的执行单元不仅能任意地读写内存，同时还能访问由软件管理的缓存，也称为共享内存。

CUDA架构的所有这些功能都是为了使GPU不仅能执行传统的图形计算，还能高效的执行通用计算。

1.4.2 CUDA架构的使用

NVIDIA专门开发了一款编译器来编译CUDA C语言。现在，用户不再需要了解OpenGL或者DirectX图形编程结构，也不需要将通用计算问题伪装为图形计算问题。

1.5 CUDA的应用

医学影像、流体力学等等

第二张 入门

2.1 本章目标

配置环境

2.2 开发环境

  • 支持CUDA的图形处理器
  • NVIDIA设备驱动程序
  • CUDA开发工具箱
  • 标准C编译器

由于CUDA C应用程序将在两个不同的处理器上执行计算，因此需要两个编译器。其中一个编译器为GPU编译代码，另一个为CPU编译代码。

第三章 CUDA简介

3.1 本章目标

第一段cuda c代码、了解host和device之间的区别、了解其他信息

3.2 第一个程序

3.2.1 Hello，World！

将CPU以及系统的内存成为主机（host）。而将GPU及其内存成为设备（device）。

在GPU设备上执行的函数通常称为核函数（kernel）。

没有核函数，只考虑在主机运行的CUDA代码和标准的C在很大程度上是没有区别的。

1
2
3
4
5

int main(void){
  printf("Hello, World!\n");
  return 0;
}

3.2.2 核函数调用

在上面的示例中添加核函数

1
2
3
4
5
6
7
8
9
10
11
#include <iostream>

__global__ void kernel(void){

}

int main(void){
  kernel<<<1,1>>>();
  printf("Hello,World!\n");
  return 0;
}

上面的代码有两处新增：

  • 一个空的核函数kernel()，并且带有修饰符__global__
  • 对这个空的核函数的调用语句，并且带有修饰字符<<<1,1>>>

cuda的host代码默认是由系统的标准C编译器来编译的（如Linux的GNU gcc和Windodws的VS），NVIDIA工具只是将代码交给host编译器，它表现出的行为就好像CUDA不存在一样。

而当遇到具有__global__修饰符的函数时，编译器就会将该函数编译为在device上运行。在此例子中，函数kernel()将被交给编译device代码的编译器，而main()函数将被交给host编译器。

而对kernel()函数的调用语句则使用了一种尖括号和两个数值的方式。这将在后面相似介绍。

3.2.3 传递参数

以下代码展示了如何像核函数传递参数并取得返回结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include "stdio.h"

__global__ void add(int a, int b, int *c){
  *c = a+b;
}

int main(){
  int c;
  int *dev_c;
  HANDLE_ERROR(cudaMalloc((void**)&dev_c, size(int)));

  add<<<1,1>>>(2,7,dev_c);

  HANDLE_ERROR(cudaMemcpy(&c,
    dev_c,
    sizeof(int),
    cudaMemcpyDeviceToHost));
   printf("2+7 = %d\n",c);
   cudaFree(dev_c);
   return 0;
}

以上多行代码包含两个概念：

  • 可以像调用C函数那样将参数传递给核函数
  • 当设备执行任何有用的操作时，都需要分配内存，例如将计算值返回给主机

cudaMalloc()函数：（注意，分配内存的指针不是该函数的返回值，这点与malloc()不同）

  • 参数一： 一个指针，指向用于保存新分配内存地址的变量。注意，由于C语言中，指针传递是本身也是值传递的，所以为了使指针本身的值（不是指针地址指向的值）可以改变，因此在传递时要使用双重指针void**，这样做的主要原因还是因为分配内存的指针最终不是通过函数返回，而是直接改变参数值导致的（如果传的是一重指针，则改变的是pd指向的内存空间的数据，而不是pd本身，所以pd也就不能指向GPU的内存了）。
  • 参数二：分配内存的大小

CUDA C的简单行及其强大功能在很大成都上都是来源于它淡化了主机代码和设备代码之间的差异。然而，程序员一定不能在主机代码中对cudaMalloc()返回的指针进行解引用（Dereference）。主机代码可以将这个指针作为参数传递，对其执行算术运算，甚至可以将其转换为另一种不同的类型。但是，绝对不饿昆虫使用这个指针来读取或写入内存。

CUDA中对设备指针的使用限制总结如下：

  • 可以将cudaMalloc()分配的指针传递给在设备上执行的函数
  • 可以在设备代码中使用cudaMalloc()分配的指针进行内存读写操作
  • 可以将cudaMalloc()分配的指针传递给在主机上执行的函数
  • 不能在主机代码中使用cudaMalloc()分配的指针进行内存读写操作

在主机代码中，可以通过调用cudaMemcpy()来访问设备上的内存。这个函数调用的行为类型与标准C中的memcpy()，只不过多了一个参数来指定设备内存指针究竟是源指针还是目标指针。如，当最后一个参数为cudaMemcpyDeviceToHost时，代表运行时源指针是一个设备指针，而目标指针是以个主机指针。此外还有参数cudaMemcpyHostToDevice和cudaMemcpyDeviceToDevice等，如果源指针和目标指针都是位于主机上，那么可以直接调用标准C的memcpy()函数。

3.3查询设备

对于拥有多个支持CUDA的设备，需要通过某种方式来确定使用的是哪一个设备。

首先，我们希望知道在系统中有多少个设备是支持CUDA架构的，并且这些设备能够运行基于CUDA C编写的核函数。要获得CUDA设备的数量，可以调用cudaGetDeviceCount()。

1
2
int count;
HANDLE_ERROR(cudaGetDeviceCount(&count));

在调用cudaGetDeviceCount()后，可以对每个设备进行迭代，并查询各个设备的相关信息。CUDA runtime将返回一个cudaDeviceProp类型的结构，其中包含了设备的相关属性。相关属性的含义可见书p20。可以利用cudaGetDeviceProperties()来获得i号设备的属性:

1
2
3
4
5
6
7
8
9
10
11
12
13

#include <iostream>

int main(){
    cudaDeviceProp prop;
    int count;
    cudaGetDeviceCount(&count);
    for(int i =0 ; i< count; i++){
	cudaGetDeviceProperties(&prop, i);
  //对设备的属性执行某些操作
    }
    std::cout<<count<<std::endl;
}

在知道了每个可用的属性以后，接下来就可以进行一些具体的操作，如：

1
std::cout<<prop.major<<std::endl;

3.4 设备属性的使用

根据在cudaGetDeviceCount()和cudaGetDeviceProperties()中返回的结果，我们可以对每个设备进行迭代，来找到我们期望的某些达到要求的设备。但是这种迭代操作执行起来有些繁琐，因此CUDA runtime提供了一种自动方式来执行这个迭代操作。首先，找出希望设备拥有的属性并将这些属性填充到一个cudaDeviceProp结构。

1
2
3
4
cudaDeviceProp prop;
memset(&prop, 0 , sizeof(cudaDeviceProp));
prop.major = 1;
prop.minor = 3;

之后，将该结构传递给cudaChooseDevice()，这样CUDA runtime运行时将查找是否存在某个设备满足这些条件，并返回一个设备ID，我们可以将这个设备ID传递给cudaSetDevice()。随后，所有的设备操作都将在这个设备上执行。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include <iostream>
using std::cout;
using std::endl;

int main(){
  cudaDeviceProp prop;
  memset(&prop, 0 , sizeof(cudaDeviceProp));
  prop.major = 1;
  prop.minor = 3;

  int dev;
  cudaGetDevice(&dev);

  cudaChooseDevice(&dev, &prop);
  cout<<"ID:"<<dev<<endl;
}

第四章 CUDA C并行编程

4.1 本章目标

  • 了解CUDA在实现并行性时采用的一种重要方式。
  • 用CUDAC编写第一段并行代码

4.2 CUDA并行编程

4.2.1 矢量求和运算

假设有两组数据，将这两组数据中对应的元素两两想加，并将结果保存在第三个数组中。

  1. 基于CPU的矢量求和
1
2
3
4
5
6
7
8
9
#define N 10

void add(int *a, int *b, int *c){
  int tid = 0; //这是第0个CPU，因此索引从0开始
  while(tid<N){
    c[tid] = a[tid] + b[tid];
    tid += 1; // 由于只有一个CPU，因此每次递增1
  }
}

  上面将代码特意写成方便修改为并行代码的形式，如，在双核处理器上，tid的设置可以分别为0和1，tid递增的大小可以改为2等。这就相当于在两个CPU上，一个执行奇数位想加，一个执行偶数位想加。

  1. 基于GPU的矢量求和

首先给出mian()函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#define N 10
int main(){
  int a[N],b[N],c[N];
  int *dev_a, *dev_b, *dev_c;

  //在GPU上分配内存，注意这里要知道为什么使用void**
  cudaMalloc( (void**)&dev_a, N*sizeof(int));
  cudaMalloc( (void**)&dev_a, N*sizeof(int));
  cudaMalloc( (void**)&dev_a, N*sizeof(int));

  ...//创建a，b数组并赋值

  //将数组a，b复制到GPU
  cudaMemcpy(dev_a, a, N*sizeof(int),cudaMemcpyHostToDevice);
  cudaMemcpy(dev_b, b, N*sizeof(int), cudaMemcpyHostToDevice);

  add<<<N,1>>>(dev_a, dev_b, dev_c);

  //将数组c从GPU复制到CPU
  cudaMemcpy(c, dev_c, N*sizeof(int), cudaMemcpyDeviceToHost);

  ...//显式结果

  //释放GPU上分配的内存
  cudaFree(dev_a);
  cudaFree(dev_b);
  cudaFree(dev_c);

  return 0;
}

看一下核函数的调用：

1
add<<<N,1>>>(dev_a, dev_b, dev_c);

尖括号中的两个数值将传递给runtime，作用是告诉runtime如何启动核函数：

  • 第一个参数：表示设备在执行款核函数时使用的并行线程块的数量。
  • 参数二：需要多少个线程格（Grid）（一格表示N个线程块的集合）

我们将每个并行执行环境都称为一个线程块（Block），对于此例，将有N个线程块在GPU上运行（N个运行核函数的副本）。

问题：如何在代码中知道当前正在运行的是哪一个线程块？

回答：利用变量blockIdx.x 。 这是一个内置变量，在CUDA runtime中已经预先定义了这个变量，无需在代码中声明，该变量中包含的值就是当前执行设备代码的线程块的索引。

接下来是GPU版本的add()函数：

1
2
3
4
5
6

__global__ void add(int *a, int *b, int *c){
  int tid = blockIdx.x;  //计算机该索引处的数据
  if(tid < N)
  c[tid] = a[tid] + b[tid];
}

当启动核函数时，我们将并行线程块的数量指定为N。这个并行线程块集合就称为一个“线程格（Grid）”。因此，此例表示我们想要一个一维的线程格，其中每个线程格包含N个线程块，每个线程块的blockInx.x的值都是不同的，cuda会为每个设备代码副本提供不同的blockInx.x。

需要注意的一点是：在启动线程块数组时，数组每一维（N）的最大数量不能超过65535。这是一种硬件限制，如过启动的线程块数量超过了这个限制，那么程序将运行失败。

4.2.2 一个有趣的示例

绘制Julia集的曲线

Julia集算法：通过下面的迭代等式对复平面中的点求值。如果在计算某个点时，迭代等式的计算结果是发散的，那么这个点就不属于Julia集合。

  1. 基于CPU的Julia集

先看main函数，它通过工具库创建了一个大小合适的位图图像，接着，将一个指向位图数据的指针传递给了核函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
//julia.cpp
#include<iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"

using namespace cv;


void kernel(Mat& M);
int julia(int x, int y);
//定义一个通用结构来保存复数值,r为实部，i为虚部
struct cuComplex{
    float r;
    float i;
    cuComplex(float a , float b ):r(a),i(b){}
    float magnitude2(){return r*r+i*i;}

    cuComplex operator*(const cuComplex& a){
	return cuComplex(r*a.r-i*a.i, i*a.r + r*a.i);
    }

    cuComplex operator+(const cuComplex& a){
	return cuComplex(a.r+r, a.i+i);
    }
};

int DIM = 800;
int main(){
  cv::Mat M(DIM,DIM,CV_8UC1);
  kernel(M);
	imshow("Test",M);   //窗口中显示图像
	imwrite("E:/灰度图.jpg",M);    //保存生成的图片
	waitKey(0); //等待按任意键后窗口自动关闭
	getchar();
  return 0;
}
void kernel(Mat& M){
  for (int y=0;y<M.rows;y++)//遍历每一行每一列并设置其像素值
	{
		for (int x=0;x<M.cols;x++)
		{
		    int juliaValue = julia(x,y);
		    M.at<uchar>(x,y)=155*juliaValue+100;
		}
	}
}

//判断函数，如果该点属于集合返回1，否则返回0
int julia(int x, int y){
    const float scale = 1.5;
    float jx = scale*(float)(DIM/2 - x)/(DIM / 2);
    float jy = scale*(float)(DIM/2 - y)/(DIM / 2);
    cuComplex c(-0.8, 0.154);
    cuComplex a(jx,jy);

    int i = 0;
    for(i = 0; i<200; i++){
      a = a*a+c;
      if(a.magnitude2() > 1000)
      return 0;
    }
    return 1;
}

下面是kernel核函数对将要绘制的所有点进行迭代，并在每次迭代时调用julia来判断该点是否属于Julia集（“是”则涂红色，“否”则涂黑色）。

该函数首先将像素坐标转换为复数空间的坐标，为了将复平面的原点定位到图像中心，代码将像素位置移动了MID/2，然后，为了确保图像的范围为-1.0到1.0，我们将图像的坐标缩放了DIM/2倍。在计算处复空间中的点之后，需要判断这个点是否属于Julia集。通过迭代判断（本示例迭代200次，在每次迭代完成后，都会判断结果是否超过某个阈值），如果属于集合，就返回1，否则，返回0。最后运行指令g++ julia.cpppkg-config —cflags —libs opencv-o julia生成的效果图如下：

  1. 基于GPU的Julia集
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include<iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"

using namespace cv;
__global__ void  kernel(unsigned char* ptr);
__device__ int julia(int x, int y);
const int DIM = 800;
struct cuComplex{
    float r;
    float i;
    __device__ cuComplex(float a , float b ):r(a),i(b){}
    __device__ float magnitude2(){return r*r+i*i;}

    __device__ cuComplex operator*(const cuComplex& a){
	return cuComplex(r*a.r-i*a.i, i*a.r + r*a.i);
    }

    __device__ cuComplex operator+(const cuComplex& a){
	return cuComplex(a.r+r, a.i+i);
    }
};
int main(){

    cv::Mat M(DIM,DIM,CV_8UC1);

    unsigned char bitmap[DIM][DIM];    
    unsigned char *dev_bitmap; //定义一定char二维数组，用来存储GPU传过来的结果
    cudaMalloc(&dev_bitmap, DIM*DIM);

    dim3 grid(DIM,DIM);
    kernel<<<grid,1>>>(dev_bitmap);

    cudaMemcpy(bitmap,dev_bitmap, DIM*DIM , cudaMemcpyDeviceToHost);
   	for (int y=0;y<M.rows;y++)        //遍历每一行每一列并设置其像素值
	{
		for (int x=0;x<M.cols;x++)
		{
		    M.at<uchar>(x,y)=bitmap[y][x];
		    //M.at<uchar>(x,y)=juliaValue+100;
		}
	}
    imshow("Test",M);   //窗口中显示图像
    waitKey(0);
    getchar();
    return 0;
}
__global__ void  kernel(unsigned char* ptr){
    int x = blockIdx.x;
    int y = blockIdx.y;
    int offset = x+ y*gridDim.x;
    int juliaValue = julia(x,y);
    ptr[offset] = 255*juliaValue;

}
__device__ int julia(int x, int y){
    const float scale = 1.5;
    float jx = scale*(float)(DIM/2 - x)/(DIM/2);
    float jy = scale*(float)(DIM/2 - y)/(DIM/2);
    cuComplex c(-0.8, 0.156);
    cuComplex a(jx,jy);
    int i = 0 ;
    for(i = 0; i<200; i++){
	a = a*a +c;
	if(a.magnitude2() > 1000)
	    return 0;
    }
    return 1;
}

这里需要注意的是，在程序中指定了多个并行线程块来执行函数kernel。并且，使用了一种新的类型来声明了一个二维的线程格：

1
dim3 grid(DIM, DIM);

类型dim3并不是标准C定义的类型，它可以表是一个三维数组，至于为什么不直接用二维数组，CUDA开发人员主要是为了日后的扩展，所以用三维数组来表示二维数组，数组的第三维默认为1。下面的代码将线程块grid传递给CUDA运行时：

1
kernel<<<grid,1>>>(dev_bitmap);

代码中还使用了修饰符__device__，这代表代码将在GPU而不是主机上运行，由于这些函数已声明为__device__，因此只能从其他__device__函数或者__global__函数中调用它们。

通常，我们将在GPU上启动的线程块集合称为一个线程格。从名字的含义可以看出，线程格既可以是一维的线程块集合，也可以是二维的线程块集合。核函数的每个副本都可以通过内置变量blockIdx来判断哪个线程块正在执行它。塌秧，还可以通过内置变量gridDim来获得线程格的大小。

第五章 线程协作

5.1 本章目标

  • 了解CUDA C中的线程
  • 了解不同线程之间的通信机制
  • 了解并行执行线程的同步机制

5.2 并行线程块的分解

当启动核函数时，我们会指定第一个参数的值，也就是指定多个并行副本，我们将这些兵行副本称为线程块（Block）。尖括号中的第二个参数表示CUDA运行时在每个线程块中创建的线程数量，因此，总共启动的线程数量可按下面的公式计算：

5.2.1 矢量求和：重新回顾

使用线程块中的并行线程，能够完成一些并行线程块无法完成的工作。

1.使用线程实现GPU上的矢量求和（即在一个线程块内设置多条线程）

相较于之前的矢量求和，需要修改两个地方，第一是将下面的代码#1式改成代码#2式：

1
2
3
add<<<N,1>>>(dev_a, dev_b, dev_C);

add<<<1,N>>>(dev_a, dev_b, dev_c);

第二是修改索引方式，有限现在只有一个线程块，所以不能再用blockIdx来获取索引了，而应使用线程索引，如下所示：

1
int tid = threadIdx.x;

完整的代码如下所是：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include<iostream>

const int N = 10;
/*#define N 20*/
__global__ void add(int * dev_a, int* dev_b, int* dev_c){
    int tid = threadIdx.x; //注意此处使用的是线程索引
    if(tid<N)
    dev_c[tid] = dev_a[tid] + dev_b[tid];
}

int main(){
    int a[N];
    int b[N];
    int c[N];

    for(int i=0;i<N;i++)
	a[i]=i;
    for(int i=0;i<N;i++)
	b[i]=i;

    int* dev_a;
    int* dev_b;
    int* dev_c;

    cudaMalloc(&dev_a, sizeof(int)*N);
    cudaMalloc(&dev_b, sizeof(int)*N);
    cudaMalloc(&dev_c, sizeof(int)*N);

    cudaMemcpy(dev_a,a,sizeof(int)*N, cudaMemcpyHostToDevice); //第一个参数为目的地址，第二个为原地址，第三个为空间大小
    cudaMemcpy(dev_b,b,sizeof(int)*N, cudaMemcpyHostToDevice);

    add<<<1,N>>>(dev_a,dev_b,dev_c); //线程块数为1，每块内的线程数为N

    cudaMemcpy(c, dev_c, sizeof(int)*N, cudaMemcpyDeviceToHost);

    for(auto x:c) std::cout<<x<<std::endl;
}

2.在GPU上对更长的矢量求和

之前在第四章我们提到，由于硬件原因，线程块的数量不能超过65535。 同样，对于启动核函数时每个线程块中的线程数量，也有一定的限制。具体来说，最大的线程数量不能超过设备属性结构中maxThreadsPerBlock域的值。对于很多图形处理器而言，这个限制值是每个线程块512个线程（目前GTX 980Ti为1024个）。通过下面的代码可以获得当前机器上的最大线程数：

1
2
3
4
5
6
7
8
9
10
11
12
13
#include <iostream>

int main(){
    cudaDeviceProp prop;
    int count;
    cudaGetDeviceCount(&count);
    for(int i =0 ; i< count; i++){
	cudaGetDeviceProperties(&prop, i);
	std::cout<<count<<std::endl;
	std::cout<<prop.maxThreadsPerBlock<<std::endl;
    }
}
//output：1024

为了通过并行对长度大于1024的矢量进行相加，必须将线程和线程块结合起来才能实现，因此仍然需要改动两个地方：核函数中的索引计算方法和核函数的调用方式：（使用多个线程块，并且每个线程块包含多个线程）

首先是修改索引计算方法：

1
2

int tid = threadIdx.x + blockIdx.x * blockDim.x;

在上面的赋值语句中使用了一个新的内置变量，blockDim。对于所有线程块来说，这个变量是一个常数，保存的是线程块中每一维的线程数量。（回顾第四章，在gridDim中保存了一个类似的值，即在线程格中每一维的线程块数量。但要知道，gridDim是二维的，blockDim是三维的，只是很少用的高维索引值）

另一处修改是核函数调用本身，为了保证最终启动的线程数量不少于预期量，可以通过一种常见的技术来对需要启动的线程块数量进行向上取整，如下所示：

1
add<<< (N+127)/128, 128 >>> (dev_a, dev_b, dev_c);

上面的代码当N不是128的整数倍时，会启动过多的线程，这时候，判断语句if (tid<N)就表现处作用了，它可以确保进行计算的线程的不会对越过数组边界的内存进行读取或写入：

1
2
if (tid < N)
  c[tid] = a[tid] + b[tid];

3.在GPU上对任意长度的矢量求和