矩阵乘法

欢迎访问我的博客首页。

---

1. 一维线程的矩阵乘法

---

  用 CUDA 多线程计算任意尺寸的矩阵 mat_a 与矩阵 mat_b 相乘，它们的尺寸分别是 $l\times m$ 和 $m\times 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>>>：

1. Dg：类型是 dim3，表示 grid 的 x、y、z 维度的长度。用于指定 block 的数量。
2. Db：类型是 dim3，表示 block 的 x、y、z 维度的长度。用于指定每个 block 内的线程数量。
3. 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

5. 参考

---

1. 矩阵乘法
2. 大数吃小数