【CUDA】 矩阵加法 Matrix Addition

矩阵加法

矩阵加法是线性代数中的一个基本运算，用于将两个大小相同的矩阵相加。加法的结果是一个与输入矩阵大小相同的矩阵。矩阵加法是逐元素进行的，天然适合CUDA并行计算。两个矩阵的加法如图1所示。

图1

线性化数组

存储空间原理请阅读文章最后一节——存储空间

至少有两种方法可以线性化二维数组。一种方法是将同一行中的所有元素连续存放，再将所有的行按行号连续存放。在图2中展示了这种行优先的存储方式。

图2

类似的，另一种线性化二维数组的方法是先连续存放同一列中的所有元素，再将所有的列按列号连续存放。FORTRAN编译器就是采用这种列优先的存储方式。请注意列优先存储方式相当于行优先存储方式的转置。希望之前有FORTRAN编程经验的人意识到CUDAC采用的是行优先存储，而不是列优先存储。同时，有很多为FORTRAN项目设计的C例程库也采用了列优先的存储方式以匹配FORTRAN编译器。这就导致了这些例程库的手册当中经常告诉使用者在C程序中调用这些例程库的时候要将输入矩阵进行转置，例如，基本线性代数子程序(BLAS)。

线程层次结构和非合并vs合并内存访问

CUDA将线程被组织成块。每个线程块最多可以有3个维度(x, y, z)。在矩阵加法中，我们只需要前两个维度(x, y)。

线程在warp内进行调度。一个warp是32个线程组成的。线程调度从x维度开始，然后是y和z维度。

每个warp以SIMD方式调度线程。这意味着一个warp中的所有线程同时执行相同的指令。为了利用这一特性，同一warp中的线程应访问相邻的内存位置。这就是所谓的合并内存访问。

Note：合并内存访问指的是warp中的线程应访问相邻的内存位置，而不是单线程连续访问相邻内存。

Code

Host代码用随机值初始化输入矩阵，并调用kernel执行矩阵加法。输入矩阵以线性化格式存储。

#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>

#include <cuda_runtime.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include "arrayAddKernels.cuh"
#include "helper_cuda.h"
#include "error.cuh"

const double EPSILON = 1.0e-10;


int main(void)
{
	int m, n, t_y, t_x;
	m = 1024;
	n = 1024;
	t_y = 32;
	t_x = 32;

	thrust::device_vector<float> A(m * n);
	thrust::device_vector<float> B(m * n);
	thrust::device_vector<float> C(m * n);

	thrust::fill(A.begin(), A.end(), 1.0f);
	thrust::fill(B.begin(), B.end(), 2.0f);

	dim3 threads(t_x, t_y);
	dim3 blocks((n + threads.x - 1) / threads.x, (m + threads.y - 1) / threads.y);

	/*thrust与cuda C指针互传。raw_pointer_cast()接受设备向量dvec的元素0的地址作为参数, 并且返回原始C指针raw_ptr。
	这个指针可用于调用CDUACAPI函数(如cudaMemset()函数)，
	或者作为参数传递到CUDACkernel函数中(如例子中的array_add_not_coalesced函数)。*/

	/*raw_pointer_cast
	creates a “raw” pointer from a pointer - like type, simply returning the wrapped pointer, should it exist.*/

	cudaEvent_t start, stop;
	checkCudaErrors(cudaEventCreate(&start));
	checkCudaErrors(cudaEventCreate(&stop));
	checkCudaErrors(cudaEventRecord(start));

	array_add_coalesced<float> <<<blocks, threads >>> (thrust::raw_pointer_cast(A.data()),
		thrust::raw_pointer_cast(B.data()),
		thrust::raw_pointer_cast(C.data()),
		m, n);

	checkCudaErrors(cudaEventRecord(stop));
	checkCudaErrors(cudaEventSynchronize(stop));
	float elapsed_time;
	checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
	printf("Time = %g ms.\n", elapsed_time);
	

	thrust::host_vector<float> h_C = C;

	bool success = true;
	for (int i = 0; i < m * n; i++) {
		if (h_C[i] != 3.0f) {
			printf("Error at (%d, %d): %f (Coalesced)\n", i / m, i % m, h_C[i]);
			success = false;
			break;
		}
	}
	if (success) printf("Coalesced: Success!\n");

	checkCudaErrors(cudaEventRecord(start));

	array_add_not_coalesced<float> <<<blocks, threads >>> (thrust::raw_pointer_cast(A.data()),
		thrust::raw_pointer_cast(B.data()),
		thrust::raw_pointer_cast(C.data()),
		m, n);

	checkCudaErrors(cudaEventRecord(stop));
	checkCudaErrors(cudaEventSynchronize(stop));
	checkCudaErrors(cudaEventElapsedTime(&elapsed_time, start, stop));
	printf("Time = %g ms.\n", elapsed_time);

	h_C = C;
	success = true;
	for (int i = 0; i < m * n; ++i)
		if (h_C[i] != 3.0f) {
			printf("Error at (%d, %d): %f (Not coalesced)\n", i / m, i % m, h_C[i]);
			success = false;
			break;
		}
	if (success) printf("Not coalesced: Success!\n");

	return 0;
}

Note：

helper_cuda.h 与error.cuh头文件为错误查验工具。后续会发布到Github。

去除checkCudaErrors等错误查验函数不影响程序运行。

非合并和合并内存访问的两个kernel如下所示。

#include <cuda_runtime.h>

template<typename T>
__global__ void array_add_not_coalesced(T* a, T* b, T* c, int m, int n)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	int j = blockIdx.y * blockDim.y + threadIdx.y;

	int index = i * m + j;

	if (i < n && j < m) {
		c[index] = a[index] + b[index];
	}
}

template<typename T>
__global__ void array_add_coalesced(T* a, T* b, T* c, int m, int n)
{
	int i = blockIdx.y * blockDim.y + threadIdx.y;
	int j = blockIdx.x * blockDim.x + threadIdx.x;

	int index = i * m + j;

	if (i < n && j < m) {
		c[index] = a[index] + b[index];
	}
}

两个kernel之间的唯一区别是前两行。

int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;

int i = blockIdx.y * blockDim.y + threadIdx.y;
int j = blockIdx.x * blockDim.x + threadIdx.x;

在第一个kernel中，第i个索引(行)由x维索引，第j个索引(列)由y维索引。这将导致相同warp中的线程访问不相邻的内存位置。由图3所示：

图3

在第二个kernel中，第i个索引(行)由y维索引，第j个索引(列)由x维索引。这将导致相同warp中的线程访问相邻的内存位置。由图4所示：

图4

下一行计算当前线程的线性化索引。

int index = i * m + j;

在最后几行中，kernel检查当前线程是否正在访问有效的内存位置并执行添加操作。

if(i < n && j < m)
    c[index] = a[index] + b[index];

性能分析

运行时间：

矩阵维度：1024 × 1024

线程块维度：32 × 32

根据运行时间分析，非合并和合并内存访问相差明显。

笔者采用设备：RTX3060 6GB

PMPP项目提供的分析：

kernel的性能是使用NvBench项目在多个gpu中测量的。研究的性能测量方法有:

内存带宽:每秒传输的数据量。

内存带宽利用率:占用内存带宽的百分比。

非合并内存访问

合并内存访问

存储空间

存储空间是现代计算机中处理器如何访问存储器的简化视图。存储空间通常与每个正在运行的应用程序相关。应用程序处理的数据和执行的指令存储在它的存储空间中，占据一定的存储单元。每个存储单元通常可容纳一个字节并且有一个唯一的地址。对于需要多个存储单元的变量存放在连续的多个存储单元中(float 变量需要4个字节、double 变量需要8 个字节)。处理器可通过提供起始地址(开始字节的存储单元地址)和字节数访问存储空间中的数据。

存储空间中的存储单元就像电话系统里的电话一样，都有一个唯一的电话号码。大多数现代计算机都有 4GB 大小以上的存储单元，1G是1073741824(2)。存储单元的地址从0开始编号，一直到最大值。由于每个存储单元都有唯一的地址，可以将存储空间看成是“扁平(flat)”的。因此，所有的多维数组需要“扁平化”为等价的一维数组。当C程序员使用多维索引语法存取多维数组中的数据时，编译器将这些操作转换成先用基地址指针指向数组的起始元素，再通过多维索引计算出偏移量，最终找到需要的存储单元。

事实上,C中的所有多维数组都是线性的。这是由于现代计算机使用的都是“扁平(flat)的存储空间。在数组静态分配的情况下，编译器允许程序员使用高维索引语法(如nums[i][i]，nums[i][i]为二维数组)的方式访问它们的元素。在底层，编译器将多维数组线性化成等价的一维数组并将多维索引翻译成一维偏移量。在数组动态分配的情况下，由于缺乏维度信息，当前的CUDA C编译器将这些转换工作留给了程序员。

参考文献：

1、大规模并行处理器编程实战（第2版）

2、PPMP