3.9 EXERCISES

1. 矩阵加法需要两个输入矩阵A和B，并产生一个输出矩阵C。输出矩阵C的每个元素都是输入矩阵A和B的相应元素的总和，即C[i][j] = A[i][j] + B[i][j]。为了简单起见，我们将只处理元素为单精度浮点数的平方矩阵。编写一个矩阵加法内核和主机stub函数，可以使用四个参数调用：指针到输出矩阵、指针到第一个输入矩阵、指针到第二个输入矩阵以及每个维度中的元素数量。按照以下说明操作：

#include <stdio.h>
#include <cuda_runtime.h>

// CUDA内核函数定义
__global__ void matrixAdd(float *C, const float *A, const float *B, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        int idx = row * N + col;
        C[idx] = A[idx] + B[idx];
    }
}

int main() {
    int N = 256; // 矩阵大小为 N x N
    size_t bytes = N * N * sizeof(float);

    // 分配主机内存
    float *h_A = (float*)malloc(bytes);
    float *h_B = (float*)malloc(bytes);
    float *h_C = (float*)malloc(bytes);

    // 初始化输入矩阵
    for (int i = 0; i < N * N; i++) {
        h_A[i] = static_cast<float>(i);
        h_B[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, bytes);
    cudaMalloc(&d_B, bytes);
    cudaMalloc(&d_C, bytes);

    // 复制数据从主机到设备
    cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);

    // 设置线程块和网格大小
    dim3 threadsPerBlock(16, 16);
    dim3 blocksPerGrid((N + threadsPerBlock.x - 1) / threadsPerBlock.x,
                       (N + threadsPerBlock.y - 1) / threadsPerBlock.y);

    // 执行内核函数
    matrixAdd<<<blocksPerGrid, threadsPerBlock>>>(d_C, d_A, d_B, N);

    // 等待CUDA完成并检查错误
    cudaDeviceSynchronize();
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(error));
        // 清理
        cudaFree(d_A);
        cudaFree(d_B);
        cudaFree(d_C);
        free(h_A);
        free(h_B);
        free(h_C);
        return -1;
    }

    // 将结果复制回主机
    cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);

    // 验证结果
    for (int i = 0; i < N * N; i++) {
        if (h_C[i] != h_A[i] + h_B[i]) {
            fprintf(stderr, "Result verification failed at element %d!\n", i);
            return -1;
        }
    }
    printf("Test PASSED\n");

    // 清理
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    free(h_A);
    free(h_B);
    free(h_C);

    return 0;
}
这段代码首先定义了一个matrixAdd内核，它计算两个矩阵的和，并将结果存储在输出矩阵中。然后在main函数中，它创建了两个输入矩阵和一个输出矩阵，将它们复制到设备内存中，调用内核函数，并将结果复制回主机内存。最后，它验证结果并清理分配的内存。

2. 矩阵-向量乘法接受输入矩阵B和向量C，并产生一个输出向量A。输出向量A的每个元素都是输入矩阵B和C的一行的点积，即 A[i] = ∑j B[i][j] + C[j]。为了简单起见，我们将只处理元素为单精度浮点数的平方矩阵。编写一个矩阵向量乘法内核和一个主机存根函数，该函数可以使用四个参数调用：指针到输出矩阵、指针到输入矩阵、指针到输入向量以及每个维度的元素数量。使用一个线程来计算输出矢量元素。

#include <stdio.h>
#include <cuda_runtime.h>

// CUDA内核函数定义
__global__ void matrixVectorMultiply(float *A, const float *B, const float *C, int N) {
    int row = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N) {
        float sum = 0.0f;
        for (int j = 0; j < N; ++j) {
            sum += B[row * N + j] * C[j];
        }
        A[row] = sum;
    }
}

int main() {
    int N = 256; // 矩阵大小为 N x N
    size_t sizeMatrix = N * N * sizeof(float);
    size_t sizeVector = N * sizeof(float);

    // 分配主机内存
    float *h_B = (float*)malloc(sizeMatrix);
    float *h_C = (float*)malloc(sizeVector);
    float *h_A = (float*)malloc(sizeVector);

    // 初始化输入矩阵B和向量C
    for (int i = 0; i < N * N; i++) {
        h_B[i] = static_cast<float>(i);
    }
    for (int i = 0; i < N; i++) {
        h_C[i] = static_cast<float>(i);
    }

    // 分配设备内存
    float *d_B, *d_C, *d_A;
    cudaMalloc(&d_B, sizeMatrix);
    cudaMalloc(&d_C, sizeVector);
    cudaMalloc(&d_A, sizeVector);

    // 复制数据从主机到设备
    cudaMemcpy(d_B, h_B, sizeMatrix, cudaMemcpyHostToDevice);
    cudaMemcpy(d_C, h_C, sizeVector, cudaMemcpyHostToDevice);

    // 设置线程块和网格大小
    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    // 执行内核函数
    matrixVectorMultiply<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

    // 等待CUDA完成并检查错误
    cudaDeviceSynchronize();
    cudaError_t error = cudaGetLastError();
    if (error != cudaSuccess) {
        fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(error));
        // 清理
        cudaFree(d_B);
        cudaFree(d_C);
        cudaFree(d_A);
        free(h_B);
        free(h_C);
        free(h_A);
        return -1;
    }

    // 将结果复制回主机
    cudaMemcpy(h_A, d_A, sizeVector, cudaMemcpyDeviceToHost);

    // 验证结果（这里跳过了验证步骤）

    // 清理
    cudaFree(d_B);
    cudaFree(d_C);
    cudaFree(d_A);
    free(h_B);
    free(h_C);
    free(h_A);

    return 0;
}

3. 如果CUDA设备的SM最多可以占用1536个线程和最多4个线程块。以下哪个块配置会导致SM中线程数量最多？
   
   答：C 。3*512

4. 对于向量加法，假设向量长度为2000，每个线程计算一个输出元素，线程块大小为512个线程。网格中将有多少个线程？
   
   答：2048

5. ***关于上一个问题，由于矢量长度的边界检查，您预计有多少warp发散？（存疑）
   
   答：A.1

因为即便在最后一个线程块中有一些线程不执行实际的向量加法操作，所有的线程都会执行相同的边界检查逻辑，因此没有warp发散。

6. 您需要编写一个在大小为400 x 900像素的图像上操作的kernel。您想为每个像素分配一个线程。您希望您的线程块是正方形的，并使用设备上每个块的最大线程数（您的设备具有计算能力3.0）。您将如何选择内核的网格尺寸和块尺寸？

对于具有计算能力3.0的NVIDIA GPU，每个线程块最多可以有1024个线程。由于您想要使用正方形的线程块，因此您应该选择一个线程块的尺寸，其边长是1024的平方根。1024的平方根是32，所以您可以选择一个线程块的尺寸为32x32，这样每个块就有$32\times 32=1024$个线程，这是最大线程数。
对于一个400x900像素的图像，您需要足够的块来覆盖所有的像素。通过将图像的尺寸除以线程块的尺寸，您可以得到网格的尺寸：

• 对于宽度（900像素），您需要 $\lceil \frac{900}{32}\rceil =29$ 个块（因为900不能被32整除，所以您需要向上取整）。
• 对于高度（400像素），您需要 $\lceil \frac{400}{32}\rceil =13$ 个块。
  因此，您的网格尺寸将是 29x13 个块，这确保了您有足够的线程块来为每个像素分配一个线程。需要注意的是，这将导致一些线程是多余的，因为29x32 = 928和13x32 = 416，而图像只有900x400。这意味着在宽度方向上将有28个线程和在高度方向上将有16个线程不会对应于图像中的任何像素。在kernel代码中，您需要添加适当的边界检查来避免这些线程进行越界操作。
  综上所述，您的内核配置应该如下：
• 线程块尺寸（blockDim）: (32, 32)
• 网格尺寸（gridDim）: (29, 13)

7. 关于上一个问题，您希望有多少个空闲线程?

为了确定有多少个空闲线程，我们需要计算网格中所有线程的总数并减去图像中像素的总数。

我们已经决定使用一个网格尺寸为29x13的线程块，每个线程块的尺寸为32x32。所以，网格中总共的线程数是：

$$\text{总线程数}=\text{网格宽度}\times \text{网格高度}\times \text{块宽度}\times \text{块高度}$$

$$\text{总线程数}=29\times 13\times 32\times 32$$

$$\text{总线程数}=387,072$$

图像中实际的像素总数是：

$$\text{像素总数}=\text{图像宽度}\times \text{图像高度}=400\times 900=360,000$$

因此，空闲线程的总数是：

$$\text{空闲线程数}=\text{总线程数}-\text{像素总数}$$

$$\text{空闲线程数}=387,072-360,000=27,072$$

所以，将会有27,072个空闲线程。这些线程应该在kernel代码中通过边界检查来处理，以确保它们不会执行任何对应于图像外部的操作。

8. 考虑一个假设的块，在到达block之前，有8个线程执行一段代码。线程需要以下时间（以微秒为单位）来执行这些部分：2.0、2.3、3.0、2.8、2.4、1.9、2.6和2.9，其余时间等待障碍。等待屏障的线程总执行时间的百分比是多少？
   答：在CUDA中，一个线程块中的所有线程在继续执行之前必须到达同步点（屏障）。这意味着所有线程必须等待直到块中最慢的线程完成其执行。在这个假设的情况下，最慢的线程需要3.0微秒来执行它的部分。

对于每个线程，等待时间是最慢线程的执行时间减去该线程的执行时间。我们需要计算每个线程的等待时间，然后将它们加起来得到总等待时间，最后计算这个总等待时间占所有线程执行时间总和的百分比。

让我们先计算总等待时间：

• 线程1的等待时间：$3.0-2.0=1.0$ 微秒
• 线程2的等待时间：$3.0-2.3=0.7$ 微秒
• 线程3的等待时间：$3.0-3.0=0.0$ 微秒（最慢的线程，没有等待时间）
• 线程4的等待时间：$3.0-2.8=0.2$ 微秒
• 线程5的等待时间：$3.0-2.4=0.6$ 微秒
• 线程6的等待时间：$3.0-1.9=1.1$ 微秒
• 线程7的等待时间：$3.0-2.6=0.4$ 微秒
• 线程8的等待时间：$3.0-2.9=0.1$ 微秒

总等待时间是所有线程等待时间的总和：

$$\text{总等待时间}=1.0+0.7+0.0+0.2+0.6+1.1+0.4+0.1=4.1\text{ 微秒}$$

总执行时间是所有线程执行时间的总和：

$$\text{总执行时间}=2.0+2.3+3.0+2.8+2.4+1.9+2.6+2.9=19.9\text{ 微秒}$$

等待时间占总执行时间的百分比是：

$$\text{百分比}=\left(\frac{\text{总等待时间}}{\text{总执行时间}}\right)\times 100$$

$$\text{百分比}=\left(\frac{4.1}{19.9}\right)\times 100\approx 20.6\%$$

因此，等待屏障的线程总执行时间的百分比大约是20.6%。

9. 指出每个多处理器可以进行以下哪些分配。在不可能的情况下，请指出限制因素（s）。
   
   好的，我们来看看每个多处理器可以进行的分配情况：

A. 在具有计算能力1.0的设备上，8个块，每个块128个线程

• 计算能力1.0的设备每个多处理器的最大线程数为768，8个块乘以128个线程等于1024个线程，超出了最大线程数的限制。因此，这种分配是不可能的，限制因素是多处理器的最大线程数。

B. 在具有计算能力1.2的设备上，8个块，每个块128个线程

• 计算能力1.2的设备每个多处理器的最大线程数为1024，8个块乘以128个线程等于1024个线程，没有超出限制。因此，这种分配是可能的。

C. 在具有计算能力3.0的设备上，8个块，每个块128个线程

• 计算能力3.0的设备每个多处理器的最大线程数为2048，8个块乘以128个线程等于1024个线程，没有超出限制。因此，这种分配是可能的。

D. 在具有计算能力1.0的设备上，16个块，每个块64个线程

• 计算能力1.0的设备每个多处理器的最大线程数为768，而16个块乘以64个线程等于1024个线程，超出了最大线程数的限制。此外，计算能力1.0的设备每个多处理器的最大块数为8，这里的16个块也超出了限制。因此，这种分配是不可能的，限制因素是多处理器的最大线程数和最大块数。

E. 在具有计算能力1.2的设备上，16个块，每个块64个线程

• 计算能力1.2的设备每个多处理器的最大线程数为1024，而16个块乘以64个线程等于1024个线程，没有超出线程数限制。但是，计算能力1.2的设备每个多处理器的最大块数为8，这里的16个块超出了限制。因此，这种分配是不可能的，限制因素是多处理器的最大块数。

F. 在具有计算能力3.0的设备上，16个块，每个块64个线程

• 计算能力3.0的设备每个多处理器的最大线程数为2048，而16个块乘以64个线程等于1024个线程，没有超出限制。计算能力3.0的设备每个多处理器的最大块数也足够容纳16个块。因此，这种分配是可能的。

10. 一位CUDA程序员说，如果他们启动一个每个块中只有32个线程的内核，他们可以在需要障碍同步的地方省略__syncthreads（）指令。你觉得这是个好主意吗？解释。
    答：不，这不是一个好主意。

__syncthreads()是CUDA编程中的一个同步屏障，它确保块内的所有线程在继续执行之前都达到这一点。这是为了避免竞争条件和数据不一致，尤其是当多个线程需要读写共享内存或者依赖于其他线程的计算结果时。

即使每个块中只有32个线程，也不能保证省略__syncthreads()指令不会导致问题。CUDA架构是SIMT（Single Instruction, Multiple Thread）架构，这意味着一组线程（称为一个线程束）将以锁步方式执行相同的指令。对于NVIDIA的CUDA架构，一个线程束通常包含32个线程，这也是所谓的warp的大小。

当一个块恰好等于一个warp的大小时，人们可能会认为每个warp内的线程是同时执行的，因此不需要同步。然而，这种假设是有风险的，因为：

1. 即使在一个warp内部，线程执行的顺序和时间仍然是不确定的，尤其是当有分支指令时。
2. 如果内核中的操作依赖于共享内存或全局内存的读写顺序，省略同步仍然可能导致错误的结果。
3. 在未来的硬件或软件更新中，warp的行为可能会改变，这使得省略__syncthreads()可能会导致代码的可移植性和可维护性问题。

因此，除非代码的逻辑确保了不同线程之间绝对不需要同步（例如，每个线程使用独立的内存位置，且不依赖于其他线程的计算结果），否则省略__syncthreads()是不安全的。通常，最好是遵循标准的同步做法，以确保代码的正确性和健壮性。

11. 一位学生提到，他能够通过使用32×32线程块的tile矩阵乘法代码来乘以两个1024×1024矩阵。他正在使用CUDA设备，每个块最多允许512个线程，每个SM最多允许8个块。他进一步提到，线程块中的每个线程都计算结果矩阵的一个元素。你会有什么反应，为什么？

答：
如果该学生声称他们使用32×32线程块来执行1024×1024矩阵的乘法，并且每个线程块中的每个线程都计算结果矩阵的一个元素，那么他们的配置与CUDA设备的硬件限制不符。具体来说，他们的线程块设置超过了每个块允许的最大线程数，因为32×32等于1024，而设备的限制是每个块最多512个线程。

这里的主要问题是，CUDA内核在执行时必须遵守硬件的限制。如果超出了这些限制，CUDA程序将无法成功执行，或者会得到错误的结果。由于CUDA内核的这个限制是由硬件决定的，因此程序员必须设计他们的线程块以适应这些限制。

因此，我的反应将是指出学生的配置错误，并建议他们重新设计线程块的尺寸。为了适应每个块最多512个线程的限制，他们可以选择如下配置之一：

• 使用16×32线程块，每块有512个线程。
• 使用32×16线程块，每块有512个线程。
• 使用其他任何乘积不超过512的线程维度配置。

此外，由于每个SM允许最多8个块，这意味着在任何给定时刻，一个SM可以同时处理的块的数量有限。这个信息对于理解内核如何在多个SM上调度是有用的，但它并不影响线程块大小的选择，只要线程块的大小遵守了每个块的最大线程数限制。