【cuda】三、矩阵相乘与coalescing writes（合并写操作）

Matrix Multiplication and Optimization

线程块

功能

• 并行执行：线程块是一组同时执行的线程。它们共同执行分配给它们的任务
• 资源共享：线程块内的线程可以共享数据和同步执行。通过共享内存（Shared Memory）和同步原语（如 __syncthreads()）实现的。
• 硬件映射：线程块的设计允许它们被有效地映射到GPU的物理硬件上。这种映射优化了执行效率，减少了线程切换和资源调度的开销。

结构

• 线程组成：一个线程块由一组线程组成，线程数量可以从1到几千不等，具体取决于CUDA架构的限制（例如，大多数CUDA设备支持每个线程块最多1024个线程）。
• 维度：线程块可以是一维、二维或三维的，这为不同类型的计算提供了灵活性。例如，二维的线程块适合于处理图像数据。
• 索引：线程块内的每个线程都有其唯一的索引，可以是一维、二维或三维的，这取决于线程块的维度。这些索引允许每个线程识别它在块内的位置，并据此处理数据。

索引和全局地址

那么有

线程索引：线程在其线程块内的二维索引 ****blockIdx的x和y。线程索引（threadIdx）表示一个线程在其所属线程块内的位置。在处理数组或矩阵时，线程索引可以用来计算要处理的元素的位置。

块索引：线程块在网格中的二维索引 blockIdx.x 和 blockIdx.y。块索引（blockIdx）表示一个线程块在整个网格（Grid）中的位置。用于确定线程块在整个问题空间中的位置。

线程块维度：blockDim.x 和 blockDim.y 表示线程块的维度。

这样就可以访问所有元素的位置地址，如果需要细节，请查看计算机组成原理课本。

例如，在二维数据处理中，一个线程的全局索引可以通过结合其线程索引和块索引来计算：

int xIndex = blockIdx.x * blockDim.x + threadIdx.x;
int yIndex = blockIdx.y * blockDim.y + threadIdx.y;

这里，blockDim.x 和 blockDim.y 表示线程块在x和y维度上的大小。通过这种方式，我们可以确定每个最小单元（thread）的地址，进行读取操作。

简单的例子：矩阵相加

例如，这里给出一个2 * 2 的线程块（Thread Blocks）。

单指令多数据（SIMD）模型

根据矩阵乘法的最基础定义公式，我们知道：结果中的每个元素的计算不依赖于结果中的其他元素。这就说明矩阵乘法任务可以进行并行。然而，我们总不能提前写好每个元素的计算公式，这样太复杂了。这就引入了SIMD模型，用于简化代码。

首先来看如下代码：

if (i < N && j < N) { // 如果这个索引在矩阵的边界内（即 i < N && j < N）
    int index = i + j * N; // 计算它的全局索引 i 和 j
    C[index] = A[index] + B[index];//独立地读取 A 和 B 中的元素，计算它们的和，然后将结果写入 C。
}

直观上来看，这就是一个串行编码中的顺序执行循环。但是，如果定义在并行的方法中，这样的串行代码就会被编译器自动转换成M*N条指令。也就是自动翻译成并行的模式。

此时**if** 语句并不是传统意义上的循环，而是一个并行执行的条件判断。

再深入一点

在更底层的层面，CUDA 运行时会将线程块分配给 GPU 上的流处理器（Streaming Multiprocessors, SMs）。SM内部包含多个CUDA核心，用于实际执行线程的计算。

线程块的调度：这个过程由CUDA运行时自动管理的，如果需要插手优化这环节，需要在核函数设计和块大小分配上间接干预。CUDA运行时会根据SM的数量和每个SM的资源情况（如寄存器、共享内存大小）来决定如何分配线程块。如果一个SM的资源不足以处理更多的线程块，新的线程块会被分配到其他SM。

每个 SM 可以同时执行多个线程，具体数量取决于 GPU 的架构和资源可用性。

• 线程调度：SMs 通过分时复用的方式在物理核心上调度线程的执行。这意味着每个核心在不同时间点可以执行不同的线程。
• 内存访问：当线程访问全局内存（如矩阵 A、B 和 C）时，存在潜在的延迟。为了最大化效率，CUDA 尝试合并对全局内存的访问，并利用局部性原理优化访问模式。
• 指令执行：GPU 采用 SIMD 或 SIMT（单指令多线程）的方式执行指令。在 SIMD 模式下，一个指令同时作用于多个数据；在 SIMT 模式下，每个线程虽然执行相同的指令序列，但可以在不同的数据上独立操作。

复杂一点：矩阵乘法

矩阵相乘是一个非常典型的例子，用于展示CUDA编程和线程块（Block）及线程（Thread）的使用。

利用tread，做矩阵乘法。

例如，这里给出一个2 * 2 的线程块（Thread Blocks）。

在这个核函数中，每个线程负责计算结果矩阵C中的一个元素。

__global__ void MatrixMultiply(float *A, float *B, float *C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {//遍历所有需要加法的地方 N 次
            sum += A[row * N + k] * B[k * N + col]; // 得到一个元素上的结果
        }
        C[row * N + col] = sum;
    }
}

主函数中调用上述核函数的方式如下：

int N = 1024; // 假设矩阵大小为1024x1024
dim3 threadsPerBlock(16, 16);
dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);

float *d_A, *d_B, *d_C;
// ... 在这里为 d_A, d_B 和 d_C 分配设备内存，并初始化数据 ...

MatrixMultiply<<<numBlocks, threadsPerBlock>>>(d_A, d_B, d_C, N);

// ... 处理完成后，从设备内存拷贝数据回主机内存，清理资源 ...

考虑性能

刚才的代码中，可以观察到两个for循环，这里可以进行优化。

“coalescing writes”（合并写操作）

“coalescing writes”（合并写操作）是一种优化内存访问模式的技术，它能显著提高内存带宽的利用效率。这种技术尤其对于全局内存访问非常重要，因为全局内存访问速度相比于核心计算速度要慢得多。

底层原理

1. 内存事务：当GPU的线程尝试访问全局内存时，这些访问被分组为内存事务。每个事务可以一次性读取或写入多个连续的字节。使用适当大小的数据类型以匹配内存事务的大小。
2. 内存对齐：为了有效地合并写操作，线程访问的内存地址应该是对齐的，并且连续线程访问的地址也应该是连续的。确保数据结构和数组在内存中对齐。
3. 线程访问模式：如果一个线程块中的所有线程都按照一定的模式（例如，线程i访问地址i）访问连续的内存地址，则这些访问可以被合并成一个或几个内存事务。设计线程块和线程索引以便线程以线性和连续的顺序访问内存。减少线程内的条件分支，以保持连续的内存访问模式。

代码

__global__ void MatrixMultiplyCoalesced(float *A, float *B, float *C, int N) {
    // 计算行和列索引
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < N; k++) {
            // 累加计算矩阵C中(row, col)位置的值
            sum += A[row * N + k] * B[k * N + col];
        }
        // 写入计算结果到矩阵C中，利用合并写操作优化
        // 每个线程按照顺序写入连续的内存地址
        C[row * N + col] = sum;
    }
}

优化点：

• 合并写操作：在写入结果到矩阵C时，每个线程写入的是连续的内存位置（C[row * N + col]）。这样，当多个线程同时写入时，由于它们访问的是连续的内存地址，这些写操作可以被合并成较少的内存事务。这种访问模式对于全局内存来说是高效的。
• 线程索引的布局：通过合理的线程索引布局（即row和col的计算方式），我们确保了线程以线性和有序的方式访问全局内存，这对于实现高效的合并写操作至关重要。