第五章 线程并行

Posted 2020-10-10 爨爨爨好

▶ 本章介绍了线程并行，并给出四个例子。长向量加法、波纹效果、点积和显示位图。

● 长向量加法（线程块并行 + 线程并行）

■ 有三个地方和上一章的单线程块并行不同，分别是 tid = threadIdx.x + blockIdx.x * blockDim.x; ； tid += blockDim.x * gridDim.x; ；以及 add <<< 128, 128 >>> (dev_a, dev_b, dev_c); 。

■ 同时使用线程块并行和线程并行，一次访问的下标范围是 gridDim.x(线程块范围) * blockDim.x(线程范围)，因此使用 tid += blockDim.x * gridDim.x; 跳到下一次访问的对应位置上去。

 1 #include <stdio.h>
 2 #include "cuda_runtime.h"
 3 #include "device_launch_parameters.h"
 4 #include "D:\\Code\\CUDA\\book\\common\\book.h"
 5 
 6 
 7 #define N   (33 * 1024)
 8 
 9 __global__ void add(int *a, int *b, int *c)
10 {
11     int tid = threadIdx.x + blockIdx.x * blockDim.x;// 与单线程块并行不同
12     while (tid < N)
13     {
14         c[tid] = a[tid] + b[tid];
15         tid += blockDim.x * gridDim.x;// 与单线程块并行不同
16     }
17     return;
18 }
19 
20 int main(void)
21 {
22     int *a, *b, *c;
23     int *dev_a, *dev_b, *dev_c;
24 
25     // 申请内存和显存
26     a = (int*)malloc(N * sizeof(int));
27     b = (int*)malloc(N * sizeof(int));
28     c = (int*)malloc(N * sizeof(int));
29     cudaMalloc((void**)&dev_a, N * sizeof(int));
30     cudaMalloc((void**)&dev_b, N * sizeof(int));
31     cudaMalloc((void**)&dev_c, N * sizeof(int));
32     
33     // 数组填充
34     for (int i = 0; i < N; i++)
35     {
36         a[i] = i;
37         b[i] = 2 * i;
38     }
39 
40     // 将内存中的a和b拷贝给显存中的dev_a和dev_b
41     cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
42     cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice);
43 
44     // 调用核函数
45     add <<< 128, 128 >>> (dev_a, dev_b, dev_c);// 与单线程块并行不同
46 
47     // 将显存中的dev_c从显存拷贝回内存中的c
48     cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost);
49 
50     // 检验结果
51     bool success = true;
52     for (int i = 0; i < N; i++)
53     {
54         if ((a[i] + b[i]) != c[i])
55         {
56             printf("Error at i==%d:\\n\\t%d + %d != %d\\n", i, a[i], b[i], c[i]);
57             success = false;
58             break;
59         }
60     }
61     if (success)
62         printf("We did it!\\n");
63 
64     // 释放内存和显存
65     free(a);
66     free(b);
67     free(c);
68     cudaFree(dev_a);
69     cudaFree(dev_b);
70     cudaFree(dev_c);
71 
72     getchar();
73     return 0;
74 }

 

● 波纹效果

■ 二维的坐标映射，将blockId.x，threadIdx.x，blockId.y，threadIdx.y映射到相应的下标上去，经常用得到。

■ 大部分技术封装到了bitmap.anim_and_exit()（接受两个函数指针，生成动画和清理显存），没有太多值得讨论的内容。

 1 #include <stdio.h>
 2 #include "cuda_runtime.h"
 3 #include "device_launch_parameters.h"
 4 #include "D:\\Code\\CUDA\\book\\common\\book.h"
 5 #include "D:\\Code\\CUDA\\book\\common\\cpu_anim.h"
 6 
 7 #define DIM 1024
 8 #define PI 3.1415926535897932f
 9 
10 struct DataBlock {
11     unsigned char   *dev_bitmap;
12     CPUAnimBitmap  *bitmap;
13 };
14 
15 __global__ void kernel(unsigned char *ptr, int ticks)//计算帧图像中每一点的灰度值
16 {
17     //标准的坐标映射
18     int x = threadIdx.x + blockIdx.x * blockDim.x;
19     int y = threadIdx.y + blockIdx.y * blockDim.y;
20     int offset = x + y * blockDim.x * gridDim.x;
21 
22     float fx = x - DIM / 2;
23     float fy = y - DIM / 2;
24     float d = sqrtf(fx * fx + fy * fy);
25     unsigned char grey = (unsigned char)(128.0f + 127.0f *cos(d / 10.0f - ticks / 7.0f) / (d / 10.0f + 1.0f));
26     ptr[offset * 4 + 0] = grey;
27     ptr[offset * 4 + 1] = grey;
28     ptr[offset * 4 + 2] = grey;
29     ptr[offset * 4 + 3] = 255;
30     
31     return;
32 }
33 
34 void generate_frame(DataBlock *d, int ticks)//生成一帧图像
35 {
36     dim3    blocks(DIM / 16, DIM / 16);
37     dim3    threads(16, 16);
38     kernel << <blocks, threads >> >(d->dev_bitmap, ticks);
39 
40     cudaMemcpy(d->bitmap->get_ptr(), d->dev_bitmap, d->bitmap->image_size(), cudaMemcpyDeviceToHost));
41 
42     return;
43 }
44 
45 void cleanup(DataBlock *d)//释放显存
46 {
47     cudaFree(d->dev_bitmap);
48 }
49 
50 int main(void)
51 {
52     DataBlock data;
53     CPUAnimBitmap bitmap(DIM, DIM, &data);
54     data.bitmap = &bitmap;
55 
56     cudaMalloc((void**)&data.dev_bitmap, bitmap.image_size());
57 
58     bitmap.anim_and_exit((void(*)(void*, int))generate_frame, (void(*)(void*))cleanup);
59 
60     getchar();
61     return 0;
62 }

■ 程序输出，动态效果，从中间向四周扩散的波动。

 

● 点积（使用共享内存）

■ 在考虑线程块大小的时候经常用到向上取整，这里使用了技巧 ceil( a / b ) == floor( (a-1) / b) + 1

■ 算法总体想法是在GPU中将很长的向量分段放入GPU的各线程块中，每个线程块利用共享内存和多线程分别计算乘法和加法。结果整理为每个线程块输出一个浮点数，置于全局内存中，这样就将待计算的元素数量降到了 gridDim.x 的水平，再返回CPU中完成剩下的加法。

■ 算法预先规定了每个线程块使用256个线程（blockDim.x == 256），那么使用的线程块数量应该满足 gridDim.x * blockDim.x ≥ N（待计算的向量长度），另外代码中规定线程块数量至少为32（？书中说“选择其他的只可能产生更高或更差的性能，这取决于CPU和GPU的相对速度”）

■ 在核函数中使用了既定大小的共享内存 __shared__ float cache[threadsPerBlock]; ，并采用 __syncthreads(); 函数进行线程同步（因为接下来要进行规约运算，前提就是该线程块内所有的线程已经独立计算完毕）。

 1 #include <stdio.h>
 2 #include "cuda_runtime.h"
 3 #include "device_launch_parameters.h"
 4 #include "D:\\Code\\CUDA\\book\\common\\book.h"
 5 
 6 #define imin(a,b) (a<b?a:b)
 7 #define sum_squares(x)  (x*(x+1)*(2*x+1)/6)//平方和计算式
 8 
 9 const int N = 33 * 1024;
10 const int threadsPerBlock = 256;
11 const int blocksPerGrid = imin(32, (N + threadsPerBlock - 1) / threadsPerBlock);
12 
13 __global__ void dot(float *a, float *b, float *c)
14 {
15     __shared__ float cache[threadsPerBlock];
16     int tid = threadIdx.x + blockIdx.x * blockDim.x;
17     int cacheIndex = threadIdx.x;
18 
19     float   temp = 0.0f;
20     while (tid < N)
21     {
22         temp += a[tid] * b[tid];
23         tid += blockDim.x * gridDim.x;
24     }
25 
26     cache[cacheIndex] = temp;//局地内存转入共享内存
27 
28     __syncthreads();//线程同步
29 
30     int i = blockDim.x / 2;//二分规约，要求每个线程块的线程数必须是2^k形式
31     while (i != 0)
32     {
33         if (cacheIndex < i)
34             cache[cacheIndex] += cache[cacheIndex + i];
35         __syncthreads();
36         i /= 2;
37     }
38 
39     if (cacheIndex == 0)//每个线程块的0号线程将，将计算结果从共享内存转入全局内存
40         c[blockIdx.x] = cache[0];
41 
42     return;
43 }
44 
45 
46 int main(void)
47 {
48     int i;
49     float   *a, *b, c, *partial_c;
50     float   *dev_a, *dev_b, *dev_partial_c;
51 
52     a = (float*)malloc(N * sizeof(float));
53     b = (float*)malloc(N * sizeof(float));
54     partial_c = (float*)malloc(blocksPerGrid * sizeof(float));
55     cudaMalloc((void**)&dev_a, N * sizeof(float));
56     cudaMalloc((void**)&dev_b, N * sizeof(float));
57     cudaMalloc((void**)&dev_partial_c, blocksPerGrid * sizeof(float));
58 
59     for (i = 0; i < N; i++)
60     {
61         a[i] = i;
62         b[i] = 2 * i;
63     }
64 
65     cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
66     cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice); 
67 
68     dot <<< blocksPerGrid, threadsPerBlock >>> (dev_a, dev_b, dev_partial_c);
69 
70     cudaMemcpy(partial_c, dev_partial_c,blocksPerGrid * sizeof(float),cudaMemcpyDeviceToHost);
71 
72     //结果在CPU中汇总
73     for (i = 0, c = 0.0f; i < blocksPerGrid; c += partial_c[i], i++);
74     printf("\\n\\tAnswer:\\t\\t%.6g\\n\\tGPU value:\\t%.6g\\n", 2 * sum_squares((float)(N - 1)), c);
75 
76     free(a);
77     free(b);
78     free(partial_c);
79     cudaFree(dev_a);
80     cudaFree(dev_b);
81     cudaFree(dev_partial_c);
82 
83     getchar();
84     return;
85 }

■ 错误的优化，想法是“只等待那些需要写入的线程来进行同步”，但是会导致有的线程无法抵达 __syncthreads() 函数而使程序停止响应。

1 while (i != 0)
2 {
3     if (cacheIndex < i)
4     {
5         cache[cacheIndex] += cache[cacheIndex + i];
6         __syncthreads();
7     }
8     i /= 2;
9 }

 ■ 正确同步的输出（左图）与不正确同步的输出（右图），共享内存中是否同步对程序结果的影响

 ■  有趣的改动，将核函数染色部分的代码改为 ptr[offset * 4 + 1] = shared[threadIdx.x][threadIdx.y]; （其他部分都不变），得到如下左图的圆形图案。特别的，如果只改 threadId.x 不改 15-threadIdx.y，得到水平方向上渐变，竖直方向上离散的右图效果。