01-CUDA编程模型与内存层次

// vector_add.cu - 向量加法示例

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

// CUDA 错误检查宏
#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            fprintf(stderr, "CUDA Error at %s:%d - %s\n", \
                    __FILE__, __LINE__, cudaGetErrorString(err)); \
            exit(EXIT_FAILURE); \
        } \
    } while(0)

// Kernel 函数：每个线程处理一个元素
__global__ void vector_add_kernel(
    const float* __restrict__ a,
    const float* __restrict__ b,
    float* __restrict__ c,
    int n
) {
    // 计算全局线程索引
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    // 边界检查
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}

// 优化版本：每个线程处理多个元素 (Grid-Stride Loop)
__global__ void vector_add_kernel_v2(
    const float* __restrict__ a,
    const float* __restrict__ b,
    float* __restrict__ c,
    int n
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;

    // Grid-Stride Loop：处理比线程数更多的元素
    for (int i = idx; i < n; i += stride) {
        c[i] = a[i] + b[i];
    }
}

// 向量化访问版本：使用 float4 一次加载 4 个元素
__global__ void vector_add_kernel_v3(
    const float4* __restrict__ a,
    const float4* __restrict__ b,
    float4* __restrict__ c,
    int n  // n 是 float4 的数量
) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;

    for (int i = idx; i < n; i += stride) {
        float4 a_val = a[i];
        float4 b_val = b[i];

        float4 c_val;
        c_val.x = a_val.x + b_val.x;
        c_val.y = a_val.y + b_val.y;
        c_val.z = a_val.z + b_val.z;
        c_val.w = a_val.w + b_val.w;

        c[i] = c_val;
    }
}

void vector_add(const float* a, const float* b, float* c, int n) {
    float *d_a, *d_b, *d_c;
    size_t size = n * sizeof(float);

    // 分配设备内存
    CUDA_CHECK(cudaMalloc(&d_a, size));
    CUDA_CHECK(cudaMalloc(&d_b, size));
    CUDA_CHECK(cudaMalloc(&d_c, size));

    // 拷贝数据到设备
    CUDA_CHECK(cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_b, b, size, cudaMemcpyHostToDevice));

    // 配置执行参数
    int block_size = 256;
    int grid_size = (n + block_size - 1) / block_size;

    // 启动 Kernel
    vector_add_kernel<<<grid_size, block_size>>>(d_a, d_b, d_c, n);
    CUDA_CHECK(cudaGetLastError());

    // 拷贝结果回主机
    CUDA_CHECK(cudaMemcpy(c, d_c, size, cudaMemcpyDeviceToHost));

    // 释放设备内存
    CUDA_CHECK(cudaFree(d_a));
    CUDA_CHECK(cudaFree(d_b));
    CUDA_CHECK(cudaFree(d_c));
}

// 1D Grid, 1D Block
__global__ void kernel_1d() {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
}

// 2D Grid, 2D Block
__global__ void kernel_2d() {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;

    // 线性索引 (row-major)
    int idx = y * (gridDim.x * blockDim.x) + x;
}

// 3D Grid, 3D Block
__global__ void kernel_3d() {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    int z = blockIdx.z * blockDim.z + threadIdx.z;

    int width = gridDim.x * blockDim.x;
    int height = gridDim.y * blockDim.y;

    int idx = z * (width * height) + y * width + x;
}

// 矩阵操作常用索引
__global__ void matrix_kernel(float* mat, int rows, int cols) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < rows && col < cols) {
        int idx = row * cols + col;  // row-major
        mat[idx] = ...;
    }
}

// 矩阵乘法 - 使用共享内存优化

#define TILE_SIZE 32

__global__ void matmul_shared(
    const float* __restrict__ A,  // [M, K]
    const float* __restrict__ B,  // [K, N]
    float* __restrict__ C,        // [M, N]
    int M, int N, int K
) {
    // 声明共享内存
    __shared__ float As[TILE_SIZE][TILE_SIZE];
    __shared__ float Bs[TILE_SIZE][TILE_SIZE];

    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    // 计算该线程负责的 C 矩阵位置
    int row = by * TILE_SIZE + ty;
    int col = bx * TILE_SIZE + tx;

    float sum = 0.0f;

    // 分块计算
    for (int t = 0; t < (K + TILE_SIZE - 1) / TILE_SIZE; t++) {
        // 协作加载 A 的一个 tile 到共享内存
        if (row < M && t * TILE_SIZE + tx < K) {
            As[ty][tx] = A[row * K + t * TILE_SIZE + tx];
        } else {
            As[ty][tx] = 0.0f;
        }

        // 协作加载 B 的一个 tile 到共享内存
        if (t * TILE_SIZE + ty < K && col < N) {
            Bs[ty][tx] = B[(t * TILE_SIZE + ty) * N + col];
        } else {
            Bs[ty][tx] = 0.0f;
        }

        // 同步，确保 tile 加载完成
        __syncthreads();

        // 计算部分和
        #pragma unroll
        for (int k = 0; k < TILE_SIZE; k++) {
            sum += As[ty][k] * Bs[k][tx];
        }

        // 同步，确保计算完成后再加载下一个 tile
        __syncthreads();
    }

    // 写回结果
    if (row < M && col < N) {
        C[row * N + col] = sum;
    }
}

/*
 * 共享内存被分成 32 个 Bank (每个 Bank 4 字节宽)
 * 同一 Warp 中的线程访问同一 Bank 的不同地址会产生 Bank Conflict
 *
 * 地址到 Bank 的映射：bank_id = (address / 4) % 32
 */

// Bank Conflict 示例
__shared__ float smem[32][32];

// 无冲突：每个线程访问不同的 Bank
float val = smem[threadIdx.x][0];  // 所有线程访问第 0 列，不同行

// 有冲突：所有线程访问同一 Bank
float val = smem[0][threadIdx.x];  // 当 threadIdx.x 是 32 的倍数时冲突

// 解决方案：Padding
__shared__ float smem_padded[32][32 + 1];  // 每行加 1 个元素的 padding

// 使用 padding 后无冲突
float val = smem_padded[0][threadIdx.x];


// 示例：转置矩阵时的 Bank Conflict 处理
__global__ void transpose_naive(float* out, const float* in, int N) {
    __shared__ float tile[32][32];

    int x = blockIdx.x * 32 + threadIdx.x;
    int y = blockIdx.y * 32 + threadIdx.y;

    // 读取 (coalesced)
    if (x < N && y < N) {
        tile[threadIdx.y][threadIdx.x] = in[y * N + x];
    }
    __syncthreads();

    // 写入 (有 bank conflict)
    x = blockIdx.y * 32 + threadIdx.x;
    y = blockIdx.x * 32 + threadIdx.y;
    if (x < N && y < N) {
        out[y * N + x] = tile[threadIdx.x][threadIdx.y];  // Bank conflict!
    }
}

__global__ void transpose_optimized(float* out, const float* in, int N) {
    __shared__ float tile[32][32 + 1];  // Padding 消除 bank conflict

    int x = blockIdx.x * 32 + threadIdx.x;
    int y = blockIdx.y * 32 + threadIdx.y;

    if (x < N && y < N) {
        tile[threadIdx.y][threadIdx.x] = in[y * N + x];
    }
    __syncthreads();

    x = blockIdx.y * 32 + threadIdx.x;
    y = blockIdx.x * 32 + threadIdx.y;
    if (x < N && y < N) {
        out[y * N + x] = tile[threadIdx.x][threadIdx.y];  // 无 bank conflict
    }
}

/*
 * 全局内存访问以 32/64/128 字节的事务为单位
 * Warp 中连续线程访问连续内存地址时可以合并为一次事务
 */

// 合并访问 (Coalesced Access) - 高效
__global__ void coalesced_access(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // 线程 0 访问 data[0]，线程 1 访问 data[1]...
        // 连续线程访问连续地址，可合并
        float val = data[idx];
    }
}

// 跨步访问 (Strided Access) - 低效
__global__ void strided_access(float* data, int n, int stride) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx * stride < n) {
        // 线程 0 访问 data[0]，线程 1 访问 data[stride]...
        // 不连续，需要多次内存事务
        float val = data[idx * stride];
    }
}

// 随机访问 - 最低效
__global__ void random_access(float* data, int* indices, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // 随机索引，几乎无法合并
        float val = data[indices[idx]];
    }
}

// 优化跨步访问：使用共享内存
__global__ void strided_access_optimized(
    float* data,
    float* output,
    int n,
    int stride
) {
    __shared__ float smem[256];

    int tid = threadIdx.x;
    int block_start = blockIdx.x * blockDim.x * stride;

    // 合并读取到共享内存
    if (block_start + tid < n) {
        smem[tid] = data[block_start + tid];
    }
    __syncthreads();

    // 从共享内存跨步读取（共享内存访问无合并要求）
    int local_idx = tid * stride % blockDim.x;
    float val = smem[local_idx];

    // 合并写回
    int out_idx = blockIdx.x * blockDim.x + tid;
    if (out_idx < n / stride) {
        output[out_idx] = val;
    }
}

#include <cuda_runtime.h>

// 1. 设备内存分配
float* d_data;
size_t size = N * sizeof(float);

cudaMalloc(&d_data, size);                    // 分配设备内存
cudaMallocManaged(&d_data, size);             // 分配统一内存
cudaMallocHost(&h_data, size);                // 分配页锁定主机内存
cudaMallocPitch(&d_data, &pitch, width, height);  // 分配 2D 内存（带对齐）

// 2. 内存拷贝
cudaMemcpy(dst, src, size, cudaMemcpyHostToDevice);
cudaMemcpy(dst, src, size, cudaMemcpyDeviceToHost);
cudaMemcpy(dst, src, size, cudaMemcpyDeviceToDevice);

// 异步拷贝
cudaMemcpyAsync(dst, src, size, kind, stream);

// 2D 拷贝
cudaMemcpy2D(dst, dpitch, src, spitch, width, height, kind);

// 3. 内存设置
cudaMemset(d_data, 0, size);                  // 按字节设置
cudaMemsetAsync(d_data, 0, size, stream);     // 异步设置

// 4. 内存释放
cudaFree(d_data);
cudaFreeHost(h_data);

// 5. 查询设备内存
size_t free_mem, total_mem;
cudaMemGetInfo(&free_mem, &total_mem);

// 统一内存示例
__global__ void kernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        data[idx] *= 2.0f;
    }
}

void unified_memory_example() {
    float* data;
    int n = 1000000;

    // 分配统一内存
    cudaMallocManaged(&data, n * sizeof(float));

    // CPU 初始化（自动在 CPU 端分页）
    for (int i = 0; i < n; i++) {
        data[i] = i;
    }

    // GPU 计算（自动迁移到 GPU）
    int block_size = 256;
    int grid_size = (n + block_size - 1) / block_size;
    kernel<<<grid_size, block_size>>>(data, n);

    // 等待完成
    cudaDeviceSynchronize();

    // CPU 访问结果（自动迁移回 CPU）
    printf("data[0] = %f\n", data[0]);

    cudaFree(data);
}

// 预取优化
void prefetch_example() {
    float* data;
    cudaMallocManaged(&data, n * sizeof(float));

    int device;
    cudaGetDevice(&device);

    // 预取到 GPU
    cudaMemPrefetchAsync(data, n * sizeof(float), device, 0);

    kernel<<<grid_size, block_size>>>(data, n);

    // 预取回 CPU
    cudaMemPrefetchAsync(data, n * sizeof(float), cudaCpuDeviceId, 0);
    cudaDeviceSynchronize();
}

// 内存建议
void memory_advise_example() {
    float* data;
    cudaMallocManaged(&data, n * sizeof(float));

    int device;
    cudaGetDevice(&device);

    // 设置为只读（优化访问）
    cudaMemAdvise(data, n * sizeof(float), cudaMemAdviseSetReadMostly, device);

    // 设置首选位置
    cudaMemAdvise(data, n * sizeof(float), cudaMemAdviseSetPreferredLocation, device);

    // 设置访问者
    cudaMemAdvise(data, n * sizeof(float), cudaMemAdviseSetAccessedBy, device);
}

// 使用页锁定内存实现 CPU-GPU 重叠
void async_transfer_example() {
    const int N = 1 << 20;
    const int num_streams = 4;
    const int chunk_size = N / num_streams;

    // 分配页锁定主机内存
    float *h_data, *h_result;
    cudaMallocHost(&h_data, N * sizeof(float));
    cudaMallocHost(&h_result, N * sizeof(float));

    // 分配设备内存
    float *d_data, *d_result;
    cudaMalloc(&d_data, N * sizeof(float));
    cudaMalloc(&d_result, N * sizeof(float));

    // 创建多个流
    cudaStream_t streams[num_streams];
    for (int i = 0; i < num_streams; i++) {
        cudaStreamCreate(&streams[i]);
    }

    // 分块异步处理
    for (int i = 0; i < num_streams; i++) {
        int offset = i * chunk_size;

        // 异步拷贝到设备
        cudaMemcpyAsync(
            d_data + offset,
            h_data + offset,
            chunk_size * sizeof(float),
            cudaMemcpyHostToDevice,
            streams[i]
        );

        // 启动 kernel
        int block_size = 256;
        int grid_size = (chunk_size + block_size - 1) / block_size;
        kernel<<<grid_size, block_size, 0, streams[i]>>>(
            d_data + offset,
            d_result + offset,
            chunk_size
        );

        // 异步拷贝回主机
        cudaMemcpyAsync(
            h_result + offset,
            d_result + offset,
            chunk_size * sizeof(float),
            cudaMemcpyDeviceToHost,
            streams[i]
        );
    }

    // 等待所有流完成
    cudaDeviceSynchronize();

    // 清理
    for (int i = 0; i < num_streams; i++) {
        cudaStreamDestroy(streams[i]);
    }
    cudaFreeHost(h_data);
    cudaFreeHost(h_result);
    cudaFree(d_data);
    cudaFree(d_result);
}

// 1. Block 内同步
__global__ void block_sync_example() {
    __shared__ float smem[256];

    // 写入共享内存
    smem[threadIdx.x] = threadIdx.x;

    // 等待所有线程完成写入
    __syncthreads();

    // 安全地读取其他线程写入的数据
    float val = smem[(threadIdx.x + 1) % 256];
}

// 2. Warp 内同步 (CUDA 9.0+)
__global__ void warp_sync_example() {
    int lane = threadIdx.x % 32;
    unsigned mask = 0xffffffff;  // 所有 32 个线程

    float val = lane;

    // Warp 内同步
    __syncwarp(mask);

    // Warp shuffle：直接从其他线程获取值
    float neighbor = __shfl_sync(mask, val, (lane + 1) % 32);

    // Warp 投票
    int vote = __ballot_sync(mask, val > 16);

    // Warp reduce
    for (int offset = 16; offset > 0; offset /= 2) {
        val += __shfl_down_sync(mask, val, offset);
    }
}

// 3. 原子操作
__global__ void atomic_example(int* counter, float* sum, float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx < n) {
        // 原子加
        atomicAdd(counter, 1);
        atomicAdd(sum, data[idx]);

        // 原子比较交换
        int old = atomicCAS(counter, 0, 1);

        // 原子最大/最小
        atomicMax(counter, idx);
        atomicMin(counter, idx);

        // 原子交换
        int prev = atomicExch(counter, idx);
    }
}

// 4. 内存栅栏
__global__ void memory_fence_example() {
    __shared__ int flag;
    __shared__ float data;

    if (threadIdx.x == 0) {
        data = 42.0f;
        __threadfence_block();  // 确保 data 对 block 内线程可见
        flag = 1;
    }

    __syncthreads();

    if (flag == 1) {
        float val = data;  // 保证读到正确的值
    }
}

#include <cooperative_groups.h>
namespace cg = cooperative_groups;

__global__ void cooperative_groups_example(float* data, int n) {
    // 获取线程块
    cg::thread_block block = cg::this_thread_block();

    // 获取 Warp
    cg::thread_block_tile<32> warp = cg::tiled_partition<32>(block);

    // Block 内同步
    block.sync();

    // Warp 内同步
    warp.sync();

    // Warp 内 reduce
    int lane = warp.thread_rank();
    float val = data[blockIdx.x * blockDim.x + threadIdx.x];

    for (int offset = warp.size() / 2; offset > 0; offset /= 2) {
        val += warp.shfl_down(val, offset);
    }

    // 只有 lane 0 有最终结果
    if (lane == 0) {
        // 使用 val
    }
}

// Grid 级别同步（需要 cooperative launch）
__global__ void grid_sync_example(float* data, int n) {
    cg::grid_group grid = cg::this_grid();

    // 第一阶段计算
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        data[idx] *= 2;
    }

    // Grid 同步：等待所有 block 完成
    grid.sync();

    // 第二阶段计算
    if (idx < n) {
        data[idx] += data[(idx + 1) % n];
    }
}

// 启动 Grid 同步 kernel
void launch_grid_sync_kernel() {
    int device;
    cudaGetDevice(&device);

    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, device);

    // 计算最大 block 数
    int num_blocks_per_sm;
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &num_blocks_per_sm, grid_sync_example, 256, 0);

    int num_blocks = num_blocks_per_sm * prop.multiProcessorCount;

    // 使用 cooperative launch
    void* args[] = { &d_data, &n };
    cudaLaunchCooperativeKernel(
        (void*)grid_sync_example,
        num_blocks, 256,
        args
    );
}

// 占用率计算
void calculate_occupancy() {
    int block_size = 256;
    int min_grid_size;
    int optimal_block_size;

    // 自动计算最优 block size
    cudaOccupancyMaxPotentialBlockSize(
        &min_grid_size,
        &optimal_block_size,
        my_kernel,
        0,  // 动态共享内存大小
        0   // block size 限制
    );

    // 计算给定 block size 的占用率
    int max_active_blocks;
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(
        &max_active_blocks,
        my_kernel,
        block_size,
        0  // 动态共享内存大小
    );

    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, 0);

    float occupancy = (float)(max_active_blocks * block_size) /
                      prop.maxThreadsPerMultiProcessor;

    printf("Occupancy: %.2f%%\n", occupancy * 100);
}

// 影响占用率的因素
/*
1. 寄存器使用量
   - 每个 SM 有固定数量的寄存器
   - 线程使用的寄存器越多，能并发的线程越少

2. 共享内存使用量
   - 每个 SM 有固定的共享内存
   - Block 使用的共享内存越多，能并发的 Block 越少

3. Block Size
   - 必须是 32 的倍数（Warp size）
   - 太小会限制并行度，太大可能受资源限制
*/

// 查看寄存器使用量
// 编译时加 --ptxas-options=-v 查看
// nvcc -ptxas-options=-v kernel.cu