欢迎访问我的博客首页。
1. 一维线程的矩阵乘法
用 CUDA 多线程计算任意尺寸的矩阵 mat_a 与矩阵 mat_b 相乘,它们的尺寸分别是 l × m l \times m l×m 和 m × n m \times n m×n。grid 和 block 的 y、z 维度长度都为 1。
1. 核函数
__global__ void matMultCUDA(float* c, const float* a, const float* b, const int m, const int n) {
const size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
const size_t row = idx / n;
const size_t col = idx % n;
float t = 0;
for (int i = 0; i < m; i++)
t += a[row * m + i] * b[i * n + col];
c[idx] = t;
}
2. 使用 C++ 调用核函数
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cmath>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
// 产生[a, b]之间的随机数。
int get_int(int a, int b) {
return rand() % (b - a + 1) + a;
}
void generate_mat(float* a, int row, int col) {
int i, j;
for (i = 0; i < row; i++)
for (j = 0; j < col; j++)
a[i * col + j] = get_int(0, 9) * 1.0;
}
void print(const char* name, const float* mat, size_t row, size_t col) {
long long sum = 0;
for (size_t i = 0; i < row; i++)
for (size_t j = 0; j < col; j++)
sum += mat[i * col + j];
printf("%s = %d.\n", name, sum);
}
int main() {
// 1.选择显卡。
int count;
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
cudaSetDevice(0);
// 2.产生输入数据。
size_t size_l = 100, size_m = 200, size_n = 300;
float* mat_a = (float*)malloc(sizeof(float) * size_l * size_m);
float* mat_b = (float*)malloc(sizeof(float) * size_m * size_n);
float* mat_c = (float*)malloc(sizeof(float) * size_l * size_n);
srand((unsigned)time(NULL));
generate_mat(mat_a, size_l, size_m);
generate_mat(mat_b, size_m, size_n);
// 3.申请显寸存放输入输出,并把数据拷贝到显存。
float* gpua, *gpub, *gpuc;
cudaMalloc((void**)&gpua, sizeof(float) * size_l * size_m);
cudaMalloc((void**)&gpub, sizeof(float) * size_m * size_n);
cudaMalloc((void**)&gpuc, sizeof(float) * size_l * size_n);
cudaMemcpy(gpua, mat_a, sizeof(float) * size_l * size_m, cudaMemcpyHostToDevice);
cudaMemcpy(gpub, mat_b, sizeof(float) * size_m * size_n, cudaMemcpyHostToDevice);
// 4.调用核函数:<<<block 数目, thread 数目, shared memory 大小>>>(参数...)。
size_t thread_total = size_l * size_n;
size_t thread_num = thread_total < 1024 ? thread_total : 1024;
size_t block_num = ceil(1.0 * thread_total / thread_num);
printf("thread_total = %d, block_num = %d, thread_num = %d.\n", thread_total, block_num, thread_num);
matMultCUDA << <block_num, thread_num, 0 >> >(gpuc, gpua, gpub, size_m, size_n);
// 5.把结果从显存复制到内存,并释放显存。
cudaMemcpy(mat_c, gpuc, sizeof(float) * size_l * size_n, cudaMemcpyDeviceToHost);
cudaFree(gpua);
cudaFree(gpub);
cudaFree(gpuc);
// 6.处理结果,释放内存。
float* mat_d = (float*)malloc(sizeof(float) * size_l * size_n);
for (size_t i = 0; i < size_l; i++)
for (size_t j = 0; j < size_n; j++) {
mat_d[i * size_n + j] = 0;
for (size_t k = 0; k < size_m; k++)
mat_d[i * size_n + j] += mat_a[i * size_m + k] * mat_b[k * size_n + j];
}
print("mat_gpu", mat_c, size_l, size_n);
print("mat_cpu", mat_d, size_l, size_n);
free(mat_a);
free(mat_b);
free(mat_c);
free(mat_d);
system("pause");
}
3. 使用 PyCUDA 调用核函数
import math
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void matMultCUDA(float* c, const float* a, const float* b, const int m, const int n) {
const size_t idx = blockIdx.x * blockDim.x + threadIdx.x;
const size_t row = idx / n;
const size_t col = idx % n;
float t = 0;
for (int i = 0; i < m; i++)
t += a[row * m + i] * b[i * n + col];
c[idx] = t;
}
""")
if __name__ == '__main__':
matMultCUDA = mod.get_function("matMultCUDA")
l, m, n = 100, 200, 300
mat_a = np.random.random(size=(l, m)).astype(np.float32)
mat_b = np.random.random(size=(m, n)).astype(np.float32)
thread_total = l * n
thread_num = thread_total if thread_total < 1024 else 1024
block_num = math.ceil(thread_total / thread_num)
dest = np.zeros(shape=(l, n), dtype=np.float32)
matMultCUDA(cuda.Out(dest), cuda.In(mat_a), cuda.In(mat_b),
np.int32(m), np.int32(n),
block=(thread_num, 1, 1), grid=(block_num, 1))
print(dest.sum())
print()
print(np.dot(mat_a, mat_b).sum())
第 22、23 行 np.random.random 产生的是 np.float64 类型的浮点数,如果不转换成 np.float32 类型,使用 C++ 的 float 类型接收这样的数会导致计算结果不正确。
第 29 行 np.int32 类型的数不能传递给 C++ size_t 类型的数。
2. 二维线程的矩阵乘法
import math
import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void matrix_mul(float *dest, float *a, float *b, int l, int m, int n) {
int i = threadIdx.x + blockDim.x * blockIdx.x;
int j = threadIdx.y + blockDim.y * blockIdx.y;
if (i >= n || j >= l)
return;
dest[j * n + i] = 0;
for (int k = 0; k < m; k++)
dest[j * n + i] += a[j * m + k] * b[k * n + i];
}
""")
if __name__ == '__main__':
matrix_mul = mod.get_function("matrix_mul")
# 1.准备数据。
l, m, n = 100, 200, 300
mat_a = np.random.uniform(low=1, high=5, size=(l, m)).astype(np.float32)
mat_b = np.random.uniform(low=1, high=5, size=(m, n)).astype(np.float32)
dest = np.zeros((l, n), dtype=np.float32)
# 2.确定grid和blick各维度的长度。
thread_total = l * n
thread_num = thread_total if thread_total < 1024 else 1024
blockDim_xy = math.ceil(math.sqrt(thread_num))
gridDim_y = math.ceil(l / blockDim_xy)
gridDim_x = math.ceil(n / blockDim_xy)
# 3.调用核函数
matrix_mul(cuda.Out(dest), cuda.In(mat_a), cuda.In(mat_b),
np.int32(l), np.int32(m), np.int32(n),
block=(blockDim_xy, blockDim_xy, 1), grid=(gridDim_x, gridDim_y))
print(dest.sum())
print()
print(np.dot(mat_a, mat_b).sum())
第 35 行 l、m、n 需要转换成 numpy 类型。第 10 行限制线程坐标,即使第 36 行创建更多线程,多余的线程也不会参加运算。
3. 大数吃小数
4. 核函数
1. 函数签名
核函数没有返回值,且定义时需要 __global__ 修饰。所以一般形式是 __global__ void Kernel(参数)。
2. 模板参数
核函数的模板参数在三层尖括号中给出 Kernel<<<Dg, Db, Ns, S>>>:
- Dg:类型是 dim3,表示 grid 的 x、y、z 维度的长度。用于指定 block 的数量。
- Db:类型是 dim3,表示 block 的 x、y、z 维度的长度。用于指定每个 block 内的线程数量。
- Ns:类型是 size_t,用于指定为每个 block 分配的共享显存大小,单位是字节。
核函数用于创建线程,所以在核函数内定义的变量是每个线程私有的。而共享显存是 1 个 block 内的所有线程共享的。
dim3 dimGrid(gridDim.x, gridDim.y);
dim3 dimBlock(blockDim.x, blockDim.y);
Kernel<<<dimGrid, dimBlock, 0>>>();
const size_t thread_x = blockIdx.x * blockDim.x + threadIdx.x; // 0 <= thread_x <= gridDim.x * blockDim.x - 1.
const size_t thread_y = blockIdx.y * blockDim.y + threadIdx.y; // 0 <= thread_y <= gridDim.y * blockDim.y - 1.
grid 和 block 是三维向量,没有指定的维度默认长度为 1。
CUDA 的 thread、block、grid 都是立体概念,也就是说它们都有三个坐标 (x,y,z)。但在本章的《矩阵相乘》与下一章的《卷积》中,我们只使用它们的两个坐标 (x,y)。这是因为我们用每个线程负责矩阵或图像平面上的某个区域,二维的线程就可以满足需求。
3. PyCUDA 调用核函数
核函数不能有 static 修饰。共享显存的大小直接在核函数内指定,共享变量不能有 extern 修饰。传递的参数类型要对应:
numpy类型 | int32 | uint64 | float32 | float64 |
---|---|---|---|---|
C++类型 | int | size_t | float | double |