▶ 本章介绍了线程并行，并给出四个例子。长向量加法、波纹效果、点积和显示位图。
● 长向量加法（线程块并行 + 线程并行）
■ 有三个地方和上一章的单线程块并行不同，分别是 tid = threadIdx.x + blockIdx.x * blockDim.x; ； tid += blockDim.x * gridDim.x; ；以及 add <<< , >>> (dev_a, dev_b, dev_c); 。
■ 同时使用线程块并行和线程并行，一次访问的下标范围是 gridDim.x(线程块范围) * blockDim.x(线程范围)，因此使用 tid += blockDim.x * gridDim.x; 跳到下一次访问的对应位置上去。
#include <stdio.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "D:\Code\CUDA\book\common\book.h" #define N (33 * 1024) __global__ void add(int *a, int *b, int *c)
{
int tid = threadIdx.x + blockIdx.x * blockDim.x;// 与单线程块并行不同
while (tid < N)
{
c[tid] = a[tid] + b[tid];
tid += blockDim.x * gridDim.x;// 与单线程块并行不同
}
return;
} int main(void)
{
int *a, *b, *c;
int *dev_a, *dev_b, *dev_c; // 申请内存和显存
a = (int*)malloc(N * sizeof(int));
b = (int*)malloc(N * sizeof(int));
c = (int*)malloc(N * sizeof(int));
cudaMalloc((void**)&dev_a, N * sizeof(int));
cudaMalloc((void**)&dev_b, N * sizeof(int));
cudaMalloc((void**)&dev_c, N * sizeof(int)); // 数组填充
for (int i = ; i < N; i++)
{
a[i] = i;
b[i] = * i;
} // 将内存中的a和b拷贝给显存中的dev_a和dev_b
cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice); // 调用核函数
add <<< , >>> (dev_a, dev_b, dev_c);// 与单线程块并行不同 // 将显存中的dev_c从显存拷贝回内存中的c
cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost); // 检验结果
bool success = true;
for (int i = ; i < N; i++)
{
if ((a[i] + b[i]) != c[i])
{
printf("Error at i==%d:\n\t%d + %d != %d\n", i, a[i], b[i], c[i]);
success = false;
break;
}
}
if (success)
printf("We did it!\n"); // 释放内存和显存
free(a);
free(b);
free(c);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c); getchar();
return ;
}
● 波纹效果
■ 二维的坐标映射，将blockId.x，threadIdx.x，blockId.y，threadIdx.y映射到相应的下标上去，经常用得到。
■ 大部分技术封装到了bitmap.anim_and_exit()（接受两个函数指针，生成动画和清理显存），没有太多值得讨论的内容。
#include <stdio.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "D:\Code\CUDA\book\common\book.h"
#include "D:\Code\CUDA\book\common\cpu_anim.h" #define DIM 1024
#define PI 3.1415926535897932f struct DataBlock {
unsigned char *dev_bitmap;
CPUAnimBitmap *bitmap;
}; __global__ void kernel(unsigned char *ptr, int ticks)//计算帧图像中每一点的灰度值
{
//标准的坐标映射
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x; float fx = x - DIM / ;
float fy = y - DIM / ;
float d = sqrtf(fx * fx + fy * fy);
unsigned char grey = (unsigned char)(128.0f + 127.0f *cos(d / 10.0f - ticks / 7.0f) / (d / 10.0f + 1.0f));
ptr[offset * + ] = grey;
ptr[offset * + ] = grey;
ptr[offset * + ] = grey;
ptr[offset * + ] = ; return;
} void generate_frame(DataBlock *d, int ticks)//生成一帧图像
{
dim3 blocks(DIM / , DIM / );
dim3 threads(, );
kernel << <blocks, threads >> >(d->dev_bitmap, ticks); cudaMemcpy(d->bitmap->get_ptr(), d->dev_bitmap, d->bitmap->image_size(), cudaMemcpyDeviceToHost)); return;
} void cleanup(DataBlock *d)//释放显存
{
cudaFree(d->dev_bitmap);
} int main(void)
{
DataBlock data;
CPUAnimBitmap bitmap(DIM, DIM, &data);
data.bitmap = &bitmap; cudaMalloc((void**)&data.dev_bitmap, bitmap.image_size()); bitmap.anim_and_exit((void(*)(void*, int))generate_frame, (void(*)(void*))cleanup); getchar();
return ;
}
■ 程序输出，动态效果，从中间向四周扩散的波动。
● 点积（使用共享内存）
■ 在考虑线程块大小的时候经常用到向上取整，这里使用了技巧 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(); 函数进行线程同步（因为接下来要进行规约运算，前提就是该线程块内所有的线程已经独立计算完毕）。
#include <stdio.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "D:\Code\CUDA\book\common\book.h" #define imin(a,b) (a<b?a:b)
#define sum_squares(x) (x*(x+1)*(2*x+1)/6)//平方和计算式 const int N = * ;
const int threadsPerBlock = ;
const int blocksPerGrid = imin(, (N + threadsPerBlock - ) / threadsPerBlock); __global__ void dot(float *a, float *b, float *c)
{
__shared__ float cache[threadsPerBlock];
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x; float temp = 0.0f;
while (tid < N)
{
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
} cache[cacheIndex] = temp;//局地内存转入共享内存 __syncthreads();//线程同步 int i = blockDim.x / ;//二分规约，要求每个线程块的线程数必须是2^k形式
while (i != )
{
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= ;
} if (cacheIndex == )//每个线程块的0号线程将，将计算结果从共享内存转入全局内存
c[blockIdx.x] = cache[]; return;
} int main(void)
{
int i;
float *a, *b, c, *partial_c;
float *dev_a, *dev_b, *dev_partial_c; a = (float*)malloc(N * sizeof(float));
b = (float*)malloc(N * sizeof(float));
partial_c = (float*)malloc(blocksPerGrid * sizeof(float));
cudaMalloc((void**)&dev_a, N * sizeof(float));
cudaMalloc((void**)&dev_b, N * sizeof(float));
cudaMalloc((void**)&dev_partial_c, blocksPerGrid * sizeof(float)); for (i = ; i < N; i++)
{
a[i] = i;
b[i] = * i;
} cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice); dot <<< blocksPerGrid, threadsPerBlock >>> (dev_a, dev_b, dev_partial_c); cudaMemcpy(partial_c, dev_partial_c,blocksPerGrid * sizeof(float),cudaMemcpyDeviceToHost); //结果在CPU中汇总
for (i = , c = 0.0f; i < blocksPerGrid; c += partial_c[i], i++);
printf("\n\tAnswer:\t\t%.6g\n\tGPU value:\t%.6g\n", * sum_squares((float)(N - )), c); free(a);
free(b);
free(partial_c);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_partial_c); getchar();
return;
}
■ 错误的优化，想法是“只等待那些需要写入的线程来进行同步”，但是会导致有的线程无法抵达 __syncthreads() 函数而使程序停止响应。
while (i != )
{
if (cacheIndex < i)
{
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
}
i /= ;
}
■ 正确同步的输出（左图）与不正确同步的输出（右图），共享内存中是否同步对程序结果的影响
■ 有趣的改动，将核函数染色部分的代码改为 ptr[offset * + ] = shared[threadIdx.x][threadIdx.y]; （其他部分都不变），得到如下左图的圆形图案。特别的，如果只改 threadId.x 不改 15-threadIdx.y，得到水平方向上渐变，竖直方向上离散的右图效果。
