国科大-高性能计算编程-作业五

作业题目

提交内容一：CUDA线程布局和内存层次

线程布局

在CUDA中，线程的布局从小到大总共有三个层次，分别为线程、线程块、网格

线程（Thread）：
- 线程是CUDA程序的最小执行单元
- 每个线程可以独立执行代码，并拥有自己的寄存器（局部内存）
线程块（Block）：
- 线程块是一定数量的线程的集合，它们可以共享一块共享内存
- 每个线程块内的线程可以通过同步原语进行协调
网格（Grid）：
- 网格是线程块的集合
- 一个CUDA内核可以启动多个线程块，所有这些块构成一个网格
- 网格可以是一维、二维或三维的结构

内存层次

CUDA在内存管理方面也存在很多个层次

寄存器（Registers）：
- 每个线程拥有自己的寄存器，这是最快的存储形式
- 用于存储局部变量和临时计算结果
共享内存（Shared Memory）：
- 由同一线程块内的所有线程共享
- 访问速度比全局内存快很多，但容量有限
- 适用于频繁访问的数据，可以显著提高性能
全局内存（Global Memory）：
- 所有线程都可以访问的内存区域
- 访问速度较慢，但容量较大
- 通常用于存储从CPU传输的数据或大规模数据集
常量内存和纹理内存（Constant and Texture Memory）：
- 常量内存用于存储不会改变的数据，对所有线程可见
- 纹理内存专用于图像处理，可以利用硬件优化的读取策略

提交内容二：使用shared memory优化代码

运行结果

原代码运行结果
引入shared memory后的运行结果

数据分析

优化后的代运行时间仅为原代码的50%

代码对比

减少全局内存访问次数：
- 原代码中的每个线程直接从全局内存读取它需要的所有数据
- 优化后的代码通过使用共享内存作为数据缓冲，每个线程块将其所需的数据加载到共享内存中，这样，线程之间可以共享这部分数据，减少了对全局内存的访问次数
合并内存访问：
- 合并访存是指多个线程同时访问位于连续地址的内存位置，这样可以最大化内存传输的效率
减少内存访问延迟：
- 共享内存比全局内存具有更低的访问延迟。通过利用共享内存缓存核心数据，每个线程处理其数据时的延迟可以显著降低，从而提高整体的执行效率。
提高带宽利用率：
- 由于减少了对全局内存的直接访问次数，并且更多的访问是合并进行的，因此提高了内存带宽的利用率

源代码

代码

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <cuda.h>
#include <cuda_runtime.h>

const int Nx = 1024 * 2;
const int Ny = 1024 * 2;
const int blockSize = 256;  // 线程块大小
const int TIMES = 10;

// 使用共享内存优化的stencil计算核函数
__global__ void two_stencil_optimized(const int n, const double * __restrict__ in_xy, double * __restrict__ out_xy) {
    extern __shared__ double tile[];
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int yindex = index / Nx;
    int xindex = index % Nx;

    // 加载输入数据到共享内存
    tile[threadIdx.x] = in_xy[index];

    __syncthreads();

    // 确保只有内部的数据才被计算
    if ((0 < xindex && xindex < (Nx - 1)) && (0 < yindex && yindex < (Ny - 1))) {
        out_xy[index] = 0.2 * (tile[threadIdx.x] + in_xy[index - 1] + in_xy[index + 1] + in_xy[index - Nx] + in_xy[index + Nx]);
    }

    __syncthreads(); // 确保所有操作在继续之前完成
}

void fill_array(const int n, double *array) {
    double init = (rand() % 1000) * 0.2;
    for (int ii = 0; ii < n; ++ii) {
        *(array + ii) = init + ii * 0.00001;
    }
}

inline int64_t GetUsec() {
    struct timeval tv;
    gettimeofday(&tv, NULL);
    return (tv.tv_sec * 1000000l) + tv.tv_usec;
}

int main() {
    srand(202405);
    double *host_in_xy = new double[Nx*Ny];
    double *host_out_xy = new double[Nx*Ny];
    fill_array(Nx*Ny, host_in_xy);
    printf("host_in_xy[1000]=%.5f\n", host_in_xy[1000]);
    
    // CUDA event 创建和计时
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    double *dev_in_xy = nullptr, *dev_out_xy = nullptr;
    cudaMalloc(&dev_in_xy, sizeof(double) * Nx * Ny);
    cudaMalloc(&dev_out_xy, sizeof(double) * Nx * Ny);
    cudaMemcpy(dev_in_xy, host_in_xy, sizeof(double) * Nx * Ny, cudaMemcpyHostToDevice);
    cudaMemset(dev_out_xy, 0, sizeof(double) * Nx * Ny);

    int numBlocks = (Nx * Ny + blockSize - 1) / blockSize;
    printf("numBlocks=%d\n", numBlocks);
    
    // warm up
    two_stencil_optimized<<<dim3(numBlocks, 1, 1), dim3(blockSize, 1, 1), blockSize * sizeof(double)>>>(Nx * Ny, dev_in_xy, dev_out_xy);
    cudaDeviceSynchronize();

    cudaEventRecord(start);
    int64_t ustart = GetUsec();
    for (int loop = 0; loop < TIMES; ++loop) {
        two_stencil_optimized<<<dim3(numBlocks, 1, 1), dim3(blockSize, 1, 1), blockSize * sizeof(double)>>>(Nx * Ny, dev_in_xy, dev_out_xy);
    }
    cudaDeviceSynchronize();
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    int64_t ufinish = GetUsec();
    cudaMemcpy(host_out_xy, dev_out_xy, sizeof(double) * Nx * Ny, cudaMemcpyDeviceToHost);
    
    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    printf("kernel time=%.5f\n", ms / TIMES);
    printf("kernel usec=%ld, host_out_xy[10000]=%.5f, host_out_xy[Nx*Ny - Nx - 16]=%.5f\n", (ufinish - ustart) / TIMES, host_out_xy[10000], host_out_xy[Nx * Ny - Nx - 16]);
    
    cudaFree(dev_in_xy);
    cudaFree(dev_out_xy);
    delete[] host_in_xy;
    delete[] host_out_xy;

    return 0;
}