欢迎访问我的博客首页。
CUDA 编程入门
1. 使用一个线程计算立方和
使用 1 个线程计算 100 个数的立方和。核函数的模板参数 block 和 grid 各个维度的长度都为 1。
1. 核函数
__global__ void sumOfSquares(int* result, const int* num, const int size) {
int sum = 0;
for (int i = 0; i < size; i++)
sum += num[i] * num[i] * num[i];
*result = sum;
}
核函数需要 global 修饰且不能有返回值，所以它有固定的头部：__global__ void。因为不能有返回值，所以只能通过参数返回，比如上面的 result。
2. 使用 C++ 调用核函数
#include <stdio.h>
#include <ctime>
#include <stdlib.h>
#include <cuda_runtime.h>
int main() {
// 1.选择显卡。
int count;
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
cudaSetDevice(0);
// 2.产生随机数。
srand((unsigned)time(NULL));
constexpr int size = 100;
int data[size];
for (int i = 0; i < size; i++)
data[i] = rand() % 10;
// 3.申请显寸存放输入输出，并把数据拷贝到显存。
int* gpudata, *result;
cudaMalloc((void**)&gpudata, sizeof(int) * size);
cudaMalloc((void**)&result, sizeof(int));
cudaMemcpy(gpudata, data, sizeof(int) * size, cudaMemcpyHostToDevice);
// 4.调用核函数：<<<block 数目, thread 数目, shared memory 大小>>>(参数...)。
sumOfSquares << <1, 1, 0 >> >(result, gpudata, size);
// 5.把结果从显存复制到内存，并释放显存。
int sum;
cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("GPUsum: %d \n", sum);
// 6.使用CPU计算的结果。
sum = 0;
for (int i = 0; i < size; i++)
sum += data[i] * data[i] * data[i];
printf("CPUsum: %d \n", sum);
system("pause");
}
第 23、24 行的 cudaMalloc 函数的第一个参数是个二重指针。该函数的作用是申请显存并用一重指针 gpudata 和 result 指向申请到的显存。注意，当实参是指针时，如果要通过函数改变实参的指向，形参须是实参的指针或引用。否则形参只是和实参指向相同地址，改变形参指向，实参不受影响。所以 cudaMalloc 函数的第一个参数是二重指针。
第 25 行从内存拷贝数据到显存。第 27 行调用核函数处理显存中的数据。第 30 行从显存拷贝数据到内存。第 31、32 行释放显存。
3. 使用 PyCUDA 调用核函数
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void sumOfSquares(int* result, const int* num, const int size) {
int sum = 0;
for (int i = 0; i < size; i++)
sum += num[i] * num[i] * num[i];
*result = sum;
}
""")
if __name__ == '__main__':
sumOfSquares = mod.get_function("sumOfSquares")
array = np.random.randint(low=1, high=5, size=[100])
dest = np.array((1), dtype=np.int32)
sumOfSquares(cuda.Out(dest), cuda.In(array),
np.int32(array.shape[0]),
block=(1, 1, 1), grid=(1, 1))
print(dest)
print()
res = 0
for i in array:
res += i * i * i
print(res)
第 20 行向核函数传递参数需要使用 numpy 类型。第 21 行指定 block 和 grid 各维度长度的参数不能使用 numpy 类型。
2. block
block 是存放线程的块，相当于一个三维空间，指定了每个维度的长度后，就指定了它容纳线程的个数。于是每个线程的 id 也可以用三维坐标表示。
2.1 使用一个 block 的多线程计算立方和
使用 32 个线程计算 1024 个数的立方和，每个线程负责 32 个数的立方和。block 的 x 维度的长度为 32，其它 2 个维度和 grid 的 3 个维度长度都为 1。
1. 核函数
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int tid = threadIdx.x;
const int size = data_size / blockDim.x;
int sum = 0;
for (int i = tid * size; i < (tid + 1) * size; i++)
sum += num[i] * num[i] * num[i];
result[tid] = sum;
}
核函数是一个线程的执行过程。比如 tid = 0 表示 0 号线程执行的是 for (int i = 0; i < size; i++)，tid = 1 表示 1 号线程执行的是 for (int i = size; i < 2 * size; i++)。sum 定义在线程中，它在 1 个线程的生命周期内有效，每个线程都有各自的 sum。而 result 是参数，被所有线程共享。
2. 使用 C++ 调用核函数
#include <ctime>
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
int main() {
// 1.选择显卡。
int count;
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
cudaSetDevice(0);
// 2.产生随机数。
srand((unsigned)time(NULL));
constexpr int size = 1024;
constexpr int thread_num = 32;
int data[size];
for (int i = 0; i < size; i++)
data[i] = rand() % 10;
// 3.申请显寸存放输入输出，并把数据拷贝到显存。
int* gpudata, *result;
cudaMalloc((void**)&gpudata, sizeof(int) * size);
cudaMalloc((void**)&result, sizeof(int) * thread_num);
cudaMemcpy(gpudata, data, sizeof(int) * size, cudaMemcpyHostToDevice);
// 4.调用核函数：<<<block 数目, thread 数目, shared memory 大小>>>(参数...)。
sumOfSquares << <1, thread_num, 0 >> >(result, gpudata, size);
// 5.把结果从显存复制到内存，并释放显存。
int sum[thread_num];
cudaMemcpy(&sum, result, sizeof(int) * thread_num, cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
int final_sum = 0;
for (int i = 0; i < thread_num; i++)
final_sum += sum[i];
printf("GPUsum: %d\n", final_sum);
// 6.使用CPU计算的结果。
final_sum = 0;
for (int i = 0; i < size; i++)
final_sum += data[i] * data[i] * data[i];
printf("CPUsum: %d \n", final_sum);
system("pause");
}
3. 使用 PyCUDA 调用核函数
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int tid = threadIdx.x;
const int size = data_size / blockDim.x;
int sum = 0;
for (int i = tid * size; i < (tid + 1) * size; i++)
sum += num[i] * num[i] * num[i];
result[tid] = sum;
}
""")
if __name__ == '__main__':
sumOfSquares = mod.get_function("sumOfSquares")
array = np.random.randint(low=1, high=3, size=[1024])
thread_num = 32
dest = np.empty(shape=(32), dtype=np.int32)
sumOfSquares(cuda.Out(dest), cuda.In(array),
np.int32(array.shape[0]),
block=(32, 1, 1), grid=(1, 1))
print(dest.sum())
print()
res = 0
for i in array:
res += i * i * i
print(res)
2.2 内存存取模式
在上面的程序中，我们使用 1 个 block 的 32 个线程处理 1024 个连续存储的数。按照核函数第 5 行的读取方式：线程 1 读取 num[0]，线程 2 读取 num[32], …, 线程 32 读取 nun[992]。然后线程 1 读取 num[1]，线程 2 读取 num[33], …, 线程 32 读取 nun[993]，直到读完。这并不是顺序读取。
顺序读取应该是：线程 1 读取 num[0]，线程 2 读取 num[1], …, 线程 32 读取 nun[31]。稍微对第 5 行的代码加以修改就能实现顺序读取。
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int tid = threadIdx.x;
int sum = 0;
for (int i = tid; i < data_size; i += blockDim.x)
sum += num[i] * num[i] * num[i];
result[tid] = sum;
}
第 4 行：线程数是 thread_num，每个线程读数据时每 thread_num 个数据读取一个，比如线程 0 读取 num[0], num[32], num[64], num[96], num[128], num[160], …, num[992]。
3. grid
grid 用于存放 block，它和 block 的关系就像 block 与线程的关系。
block(x1, y1, z1) 每个维度的长度代表这个维度上线程的个数，grid(x2, y2, z2) 每个维度的长度代表这个维度上 block 的个数。如果把线程类比为毫米单位，则 block 类似于厘米单位，grid 类似于分米单位。于是线程三个维度的长度就是 (x1x2, y1y2, z1z2)，共有 x1x2y1y2z1z2 个线程。
3.1 使用多个 block 的多线程计算立方和
上面的程序在核函数模板中指定 grid 的各维度都是 1，也就是使用 1 个 block，包含 32 个线程。一个 block 最多有 1024 个线程，当线程多时，可以使用多个 block。下面的程序中 grid 的 x 维度长度为 4，block 的 x 维度长度为 32，其它维度长度都为 1。于是 x 维度上共有 128 个线程，每个线程负责计算 8 个数的立方和。
1. 核函数
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int bid = blockIdx.x;
const int tid = threadIdx.x;
int sum = 0;
for (int i = bid * blockDim.x + tid; i < data_size; i += gridDim.x * blockDim.x)
sum += num[i] * num[i] * num[i];
result[bid * blockDim.x + tid] = sum;
}
第 2、3 行获取当前线程在 grid 的 x 维度上的坐标和在 block 的 x 维度上的坐标。
grid 的坐标可以类比为厘米单位，block 的坐标可以类比为毫米单位。blockIdx.x 是线程的厘米坐标，threadIdx.x 是线程在一个厘米单位内的毫米坐标。因为每个 block 包含 thread_num 个线程，所以线程坐标是 blockIdx.x * blockDim.x + threadIdx.x。
2. 使用 C++ 调用核函数
#include <ctime>
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
int main() {
// 1.选择显卡。
int count;
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
cudaSetDevice(0);
// 2.产生随机数。
srand((unsigned)time(NULL));
constexpr int size = 1024;
constexpr int block_num = 4;
constexpr int thread_num = 32;
int data[size];
for (int i = 0; i < size; i++)
data[i] = rand() % 10;
// 3.申请显寸存放输入输出，并把数据拷贝到显存。
int* gpudata, *result;
cudaMalloc((void**)&gpudata, sizeof(int) * size);
cudaMalloc((void**)&result, sizeof(int) * block_num * thread_num);
cudaMemcpy(gpudata, data, sizeof(int) * size, cudaMemcpyHostToDevice);
// 4.调用核函数：<<<block 数目, thread 数目, shared memory 大小>>>(参数...)。
sumOfSquares << <block_num, thread_num, 0 >> >(result, gpudata, size);
// 5.把结果从显存复制到内存，并释放显存。
int sum[block_num * thread_num];
cudaMemcpy(&sum, result, sizeof(int) * block_num * thread_num, cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
int final_sum = 0;
for (int i = 0; i < block_num * thread_num; i++)
final_sum += sum[i];
printf("GPUsum: %d\n", final_sum);
// 6.使用CPU计算的结果。
final_sum = 0;
for (int i = 0; i < size; i++)
final_sum += data[i] * data[i] * data[i];
printf("CPUsum: %d \n", final_sum);
system("pause");
}
3. 使用 PyCUDA 调用核函数
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int bid = blockIdx.x;
const int tid = threadIdx.x;
int sum = 0;
for (int i = bid * blockDim.x + tid; i < data_size; i += gridDim.x * blockDim.x)
sum += num[i] * num[i] * num[i];
result[bid * blockDim.x + tid] = sum;
}
""")
if __name__ == '__main__':
sumOfSquares = mod.get_function("sumOfSquares")
array = np.random.randint(low=1, high=3, size=[1024])
block_num = 4
thread_num = 32
dest = np.empty(shape=(block_num * thread_num), dtype=np.int32)
sumOfSquares(cuda.Out(dest), cuda.In(array),
np.int32(array.shape[0]),
block=(32, 1, 1), grid=(block_num, 1))
print(dest.sum())
print()
res = 0
for i in array:
res += i * i * i
print(res)
3.2 共享内存
上面的程序使用 block 把线程数增加到 128，然后在第 33 行把 128 个线程的计算结果拷贝到内存，在第 37、38 行处理这些结果。如果线程较多、每个线程计算的结果较大，拷贝会造成很大开销，处理这些结果也会给 CPU 带来很大负担。
因为一个 block 内的线程可以相互同步和通信，所以可以把每个 block 内的结果在 GPU 内处理后，再拷贝到内存处理。这样，需要拷贝和需要 CPU 处理的数据量就从 128 减少到 4。
1. 核函数
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int bid = blockIdx.x;
const int tid = threadIdx.x;
extern __shared__ int shared[];
shared[tid] = 0;
for (int i = bid * blockDim.x + tid; i < data_size; i += gridDim.x * blockDim.x)
shared[tid] += num[i] * num[i] * num[i];
__syncthreads();
if (tid == 0) {
for (int i = 1; i < blockDim.x; i++)
shared[0] += shared[i];
result[bid] = shared[0];
}
}
第 4 行的共享内存 shared 由 1 个 block 内的所有线程共享，其大小由调用核函数时的第 3 个模板参数指定。
2. 使用 C++ 调用核函数
#include <ctime>
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
int main() {
// 1.选择显卡。
int count;
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
cudaSetDevice(0);
// 2.产生随机数。
srand((unsigned)time(NULL));
constexpr int size = 1024;
constexpr int block_num = 4;
constexpr int thread_num = 32;
int data[size];
for (int i = 0; i < size; i++)
data[i] = rand() % 10;
// 3.申请显寸存放输入输出，并把数据拷贝到显存。
int* gpudata, *result;
cudaMalloc((void**)&gpudata, sizeof(int) * size);
cudaMalloc((void**)&result, sizeof(int) * block_num);
cudaMemcpy(gpudata, data, sizeof(int) * size, cudaMemcpyHostToDevice);
// 4.调用核函数：<<<block 数目, thread 数目, shared memory 大小>>>(参数...)。
sumOfSquares << <block_num, thread_num, thread_num * sizeof(int) >> >(result, gpudata, size);
// 5.把结果从显存复制到内存，并释放显存。
int sum[block_num];
cudaMemcpy(&sum, result, sizeof(int) * block_num, cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
int final_sum = 0;
for (int i = 0; i < block_num; i++)
final_sum += sum[i];
printf("GPUsum: %d\n", final_sum);
// 6.使用CPU计算的结果。
final_sum = 0;
for (int i = 0; i < size; i++)
final_sum += data[i] * data[i] * data[i];
printf("CPUsum: %d \n", final_sum);
system("pause");
}
第 27 行只需为每个 block 申请一个 int 大小的内存空间存放结果。
3. 使用 PyCUDA 调用核函数
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int bid = blockIdx.x;
const int tid = threadIdx.x;
__shared__ int shared[128];
shared[tid] = 0;
for (int i = bid * blockDim.x + tid; i < data_size; i += gridDim.x * blockDim.x)
shared[tid] += num[i] * num[i] * num[i];
__syncthreads();
if (tid == 0) {
for (int i = 1; i < blockDim.x; i++)
shared[0] += shared[i];
result[bid] = shared[0];
}
}
""")
if __name__ == '__main__':
sumOfSquares = mod.get_function("sumOfSquares")
array = np.random.randint(low=1, high=3, size=[1024])
block_num = 4
thread_num = 32
dest = np.empty(shape=(block_num), dtype=np.int32)
sumOfSquares(cuda.Out(dest), cuda.In(array),
np.int32(array.shape[0]),
block=(32, 1, 1), grid=(block_num, 1))
print(dest.sum())
print()
res = 0
for i in array:
res += i * i * i
print(res)
注意，使用 PyCUDA 调用核函数时，共享内存不能使用 extern 修饰。
3.3 树状加法
上面的代码通过在 block 内处理线程的结果，大大减少了需要拷贝到内存和需要 CPU 处理的数据量。但使用的方法是把一个 block 内的多个线程的结果累加到第一个线程的结果上，这个累加过程不是并行的。下面使用树状加法把这个过程并行化。
| offset | 序号0 | 序号1 | 序号2 | 序号3 | 序号4 | 序号5 | 序号6 | 序号7 | 序号8 | 序号9 | 序号10 | 序号11 | 序号12 | 序号13 | 序号14 | 序号15 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
| 2 | 1 | 1 | 5 | 3 | 9 | 5 | 13 | 7 | 17 | 9 | 21 | 11 | 25 | 13 | 29 | 15 |
| 4 | 6 | 1 | 5 | 3 | 22 | 5 | 13 | 7 | 38 | 9 | 21 | 11 | 54 | 13 | 29 | 15 |
| 8 | 28 | 1 | 5 | 3 | 22 | 5 | 13 | 7 | 92 | 9 | 21 | 11 | 54 | 13 | 29 | 15 |
| 120 |
红色的线程序号 tid 能被 2 × o f f s e t 2 \times offset 2×offset 整除，于是就把它后面偏移量为 offset 的绿色结果累加给它。比如 8 能被 2 × 4 2 \times 4 2×4 整除，就把 tid12 的结果 54 累加到 38 得到 92。
__global__ void sumOfSquares(int* result, const int* num, const int data_size) {
const int bid = blockIdx.x;
const int tid = threadIdx.x;
extern __shared__ int shared[];
shared[tid] = 0;
for (int i = bid * blockDim.x + tid; i < data_size; i += gridDim.x * blockDim.x)
shared[tid] += num[i] * num[i] * num[i];
__syncthreads();
// 树状加法。
int offset = 1;
while (offset < blockDim.x) {
if ((tid & (2 * offset - 1)) == 0) // if (tid % (2 * offset) == 0)
shared[tid] += shared[tid + offset];
offset *= 2;
__syncthreads();
}
if (tid == 0)
result[bid] = shared[0];
}
第 10 至第 16 行实现树状加法。第 12 行： x % 2 n = x & ( 2 n − 1 ) x \% 2^n = x \& (2^n - 1) x%2n=x&(2n−1)。 x % m = x & ( m − 1 ) x \% m = x \& (m - 1) x%m=x&(m−1) 的条件是 m 是 2 的整数次方。
4. 异常捕获
在 cuda 函数 cudaMalloc、cudaMemcpy 等后面使用下面的语句可以捕获异常信息。
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess)
cout << cudaGetErrorString(err) << endl;
也可以把 cuda 函数放在下面的 log 函数括号中，仅限能返回错误信息的函数：
void print(cudaError_t err, const char* file, const char* func, const int line, const char* timestamp) {
if (err != cudaSuccess)
cout << "-- file " << file << ", function " << func << ", line " << line << ", "\
<< timestamp << ":" << endl << " " << cudaGetErrorString(err) << endl;
}
#define log(msg) (print(msg, __FILE__, __FUNCTION__, __LINE__, __TIMESTAMP__))
// 使用
log(cudaMalloc(...));
4704
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
