CUDA 并行算法实战笔记 + C++示例代码（第六章重点）

#include <cuda_runtime.h>
#include <iostream>
#include <vector>
#include <chrono>
#include <cmath>

using namespace std;
using namespace chrono;

#define CHECK_CUDA_ERR(err) \
    if (err != cudaSuccess) { \
        cerr << "CUDA错误：" << cudaGetErrorString(err) << "（行号：" << __LINE__ << "）" << endl; \
        exit(1); \
    }

// 分块大小（Tile Size，适配共享内存和Warp大小）
const int TILE_SIZE = 32;
// 矩阵尺寸（必须是TILE_SIZE的整数倍）
const int MATRIX_N = 1024;

// CPU串行矩阵乘法（A[M×K] × B[K×N] = C[M×N]）
void matrixMulCPU(const float* h_A, const float* h_B, float* h_C, int M, int K, int N) {
    for (int i = 0; i < M; ++i) {
        for (int j = 0; j < N; ++j) {
            float sum = 0.0f;
            for (int k = 0; k < K; ++k) {
                sum += h_A[i * K + k] * h_B[k * N + j];
            }
            h_C[i * N + j] = sum;
        }
    }
}

// GPU核函数：共享内存分块优化矩阵乘法
__global__ void matrixMulSharedGPU(const float* d_A, const float* d_B, float* d_C, int M, int K, int N) {
    // 共享内存：存储A的Tile和B的Tile
    __shared__ float shared_A[TILE_SIZE][TILE_SIZE];
    __shared__ float shared_B[TILE_SIZE][TILE_SIZE];

    // 2D全局索引（映射到C矩阵的行和列）
    int row = threadIdx.y + blockIdx.y * blockDim.y;
    int col = threadIdx.x + blockIdx.x * blockDim.x;

    float result = 0.0f;

    // 分块遍历K维度（每次加载一个Tile到共享内存）
    for (int tile_k = 0; tile_k < (K + TILE_SIZE - 1) / TILE_SIZE; ++tile_k) {
        // 加载A的当前Tile到共享内存（避免越界）
        if (row < M && (tile_k * TILE_SIZE + threadIdx.x) < K) {
            shared_A[threadIdx.y][threadIdx.x] = d_A[row * K + tile_k * TILE_SIZE + threadIdx.x];
        } else {
            shared_A[threadIdx.y][threadIdx.x] = 0.0f; // 超出范围填0
        }

        // 加载B的当前Tile到共享内存（避免越界）
        if (col < N && (tile_k * TILE_SIZE + threadIdx.y) < K) {
            shared_B[threadIdx.y][threadIdx.x] = d_B[(tile_k * TILE_SIZE + threadIdx.y) * N + col];
        } else {
            shared_B[threadIdx.y][threadIdx.x] = 0.0f; // 超出范围填0
        }

        // 等待所有线程加载完成（块内同步）
        __syncthreads();

        // 计算当前Tile的部分和（使用共享内存，合并访问）
        for (int k = 0; k < TILE_SIZE; ++k) {
            result += shared_A[threadIdx.y][k] * shared_B[k][threadIdx.x];
        }

        // 等待所有线程计算完成，避免覆盖共享内存
        __syncthreads();
    }

    // 存储结果到C矩阵
    if (row < M && col < N) {
        d_C[row * N + col] = result;
    }
}

int main() {
    const int M = MATRIX_N; // A行数
    const int K = MATRIX_N; // A列数 = B行数
    const int N = MATRIX_N; // B列数
    const size_t size_A = M * K * sizeof(float);
    const size_t size_B = K * N * sizeof(float);
    const size_t size_C = M * N * sizeof(float);

    // 1. 主机内存初始化（A、B全1，C全0）
    vector<float> h_A(M * K, 1.0f);
    vector<float> h_B(K * N, 1.0f);
    vector<float> h_C_CPU(M * N, 0.0f);
    vector<float> h_C_GPU(M * N, 0.0f);

    // 2. CPU串行计算（基准测试）
    auto cpu_start = high_resolution_clock::now();
    matrixMulCPU(h_A.data(), h_B.data(), h_C_CPU.data(), M, K, N);
    auto cpu_time = duration_cast<seconds>(high_resolution_clock::now() - cpu_start).count();
    cout << "CPU矩阵乘法（" << MATRIX_N << "×" << MATRIX_N << "）耗时：" << cpu_time << " s" << endl;

    // 3. GPU并行计算（共享内存优化）
    float *d_A, *d_B, *d_C;
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_A, size_A));
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_B, size_B));
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_C, size_C));

    // 2D线程配置（TILE_SIZE×TILE_SIZE线程/块）
    const dim3 BLOCK_DIM(TILE_SIZE, TILE_SIZE);
    const dim3 GRID_DIM((N + BLOCK_DIM.x - 1) / BLOCK_DIM.x,
                        (M + BLOCK_DIM.y - 1) / BLOCK_DIM.y);

    // GPU计时
    cudaEvent_t gpu_start, gpu_stop;
    CHECK_CUDA_ERR(cudaEventCreate(&gpu_start));
    CHECK_CUDA_ERR(cudaEventCreate(&gpu_stop));
    CHECK_CUDA_ERR(cudaEventRecord(gpu_start, 0));

    // 数据拷贝（主机→设备）
    CHECK_CUDA_ERR(cudaMemcpy(d_A, h_A.data(), size_A, cudaMemcpyHostToDevice));
    CHECK_CUDA_ERR(cudaMemcpy(d_B, h_B.data(), size_B, cudaMemcpyHostToDevice));

    // 启动核函数（共享内存优化）
    matrixMulSharedGPU<<<GRID_DIM, BLOCK_DIM>>>(d_A, d_B, d_C, M, K, N);
    CHECK_CUDA_ERR(cudaGetLastError());

    // 数据拷贝（设备→主机）
    CHECK_CUDA_ERR(cudaMemcpy(h_C_GPU.data(), d_C, size_C, cudaMemcpyDeviceToHost));

    // 计时结束
    CHECK_CUDA_ERR(cudaEventRecord(gpu_stop, 0));
    CHECK_CUDA_ERR(cudaEventSynchronize(gpu_stop));
    float gpu_time_ms;
    CHECK_CUDA_ERR(cudaEventElapsedTime(&gpu_time_ms, gpu_start, gpu_stop));
    float gpu_time_s = gpu_time_ms / 1000.0f;

    // 4. 结果验证（前10个元素，预期值均为MATRIX_N）
    bool valid = true;
    for (int i = 0; i < 10; ++i) {
        if (abs(h_C_CPU[i] - MATRIX_N) > 1e-5 || abs(h_C_GPU[i] - MATRIX_N) > 1e-5) {
            valid = false;
            break;
        }
    }

    // 输出结果
    cout << "GPU矩阵乘法（共享内存优化）耗时：" << gpu_time_s << " s" << endl;
    cout << "结果验证：" << (valid ? "正确" : "错误") << endl;
    cout << "GPU加速比：" << (double)cpu_time / gpu_time_s << "x" << endl;

    // 5. 释放资源
    CHECK_CUDA_ERR(cudaFree(d_A));
    CHECK_CUDA_ERR(cudaFree(d_B));
    CHECK_CUDA_ERR(cudaFree(d_C));
    CHECK_CUDA_ERR(cudaEventDestroy(gpu_start));
    CHECK_CUDA_ERR(cudaEventDestroy(gpu_stop));

    return 0;
}

#include <cuda_runtime.h>
#include <iostream>
#include <vector>
#include <chrono>
#include <cmath>

using namespace std;
using namespace chrono;

#define CHECK_CUDA_ERR(err) \
    if (err != cudaSuccess) { \
        cerr << "CUDA错误：" << cudaGetErrorString(err) << "（行号：" << __LINE__ << "）" << endl; \
        exit(1); \
    }

// 积分参数：f(x) = sin(x)，积分区间[0, π]，采样点数量
const int N_SAMPLES = 10'000'000; // 1000万个采样点
const float X_START = 0.0f;
const float X_END = M_PI;
const float INTERVAL = (X_END - X_START) / N_SAMPLES; // 采样区间宽度

// CPU串行梯形法积分
float integralCPU(const vector<float>& samples) {
    float sum = 0.0f;
    for (int i = 0; i < N_SAMPLES - 1; ++i) {
        sum += 0.5f * (samples[i] + samples[i + 1]) * INTERVAL;
    }
    return sum;
}

// GPU核函数：梯形法积分 + 块内并行归约
__global__ void trapezoidalReductionGPU(const float* d_samples, float* d_partial_sums, int n) {
    extern __shared__ float shared_sum[]; // 共享内存（核函数启动时指定大小）
    int tid = threadIdx.x;
    int global_idx = tid + blockIdx.x * blockDim.x;

    // 1. 计算单个梯形面积（局部结果）
    float local_area = 0.0f;
    if (global_idx < n - 1) {
        local_area = 0.5f * (d_samples[global_idx] + d_samples[global_idx + 1]) * INTERVAL;
    }

    // 2. 局部结果存入共享内存
    shared_sum[tid] = local_area;
    __syncthreads(); // 等待块内所有线程写入完成

    // 3. 块内树状归约（步长减半，合并部分和）
    for (int stride = blockDim.x / 2; stride > 0; stride /= 2) {
        if (tid < stride) {
            shared_sum[tid] += shared_sum[tid + stride];
        }
        __syncthreads(); // 每轮合并后同步
    }

    // 4. 块内归约结果存入全局内存（线程0负责）
    if (tid = 0) {
        d_partial_sums[blockIdx.x] = shared_sum[0];
    }
}

int main() {
    // 1. 生成采样数据：f(x) = sin(x)，x ∈ [0, π]
    vector<float> h_samples(N_SAMPLES);
    for (int i = 0; i < N_SAMPLES; ++i) {
        float x = X_START + i * INTERVAL;
        h_samples[i] = sin(x);
    }

    // 2. CPU串行积分
    auto cpu_start = high_resolution_clock::now();
    float cpu_result = integralCPU(h_samples);
    auto cpu_time = duration_cast<milliseconds>(high_resolution_clock::now() - cpu_start).count();
    cout << "CPU梯形法积分结果：" << cpu_result << "，耗时：" << cpu_time << " ms" << endl;

    // 3. GPU并行积分+归约
    const int BLOCK_SIZE = 256;
    const int GRID_SIZE = (N_SAMPLES + BLOCK_SIZE - 1) / BLOCK_SIZE;
    const size_t samples_size = N_SAMPLES * sizeof(float);
    const size_t partial_sums_size = GRID_SIZE * sizeof(float);

    // 设备内存分配
    float *d_samples, *d_partial_sums;
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_samples, samples_size));
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_partial_sums, partial_sums_size));

    // 主机内存分配（存储部分和）
    vector<float> h_partial_sums(GRID_SIZE);

    // GPU计时
    cudaEvent_t gpu_start, gpu_stop;
    CHECK_CUDA_ERR(cudaEventCreate(&gpu_start));
    CHECK_CUDA_ERR(cudaEventCreate(&gpu_stop));
    CHECK_CUDA_ERR(cudaEventRecord(gpu_start, 0));

    // 数据拷贝（主机→设备）
    CHECK_CUDA_ERR(cudaMemcpy(d_samples, h_samples.data(), samples_size, cudaMemcpyHostToDevice));

    // 启动核函数（共享内存大小=BLOCK_SIZE×sizeof(float)）
    trapezoidalReductionGPU<<<GRID_SIZE, BLOCK_SIZE, BLOCK_SIZE * sizeof(float)>>>(
        d_samples, d_partial_sums, N_SAMPLES);
    CHECK_CUDA_ERR(cudaGetLastError());

    // 数据拷贝（设备→主机：部分和）
    CHECK_CUDA_ERR(cudaMemcpy(h_partial_sums.data(), d_partial_sums, partial_sums_size, cudaMemcpyDeviceToHost));

    // 4. CPU二次归约（合并块内部分和，开销极小）
    float gpu_result = 0.0f;
    for (float s : h_partial_sums) {
        gpu_result += s;
    }

    // 计时结束
    CHECK_CUDA_ERR(cudaEventRecord(gpu_stop, 0));
    CHECK_CUDA_ERR(cudaEventSynchronize(gpu_stop));
    float gpu_time;
    CHECK_CUDA_ERR(cudaEventElapsedTime(&gpu_time, gpu_start, gpu_stop));

    // 5. 结果验证（理论值=2.0）
    bool valid = abs(gpu_result - 2.0f) < 1e-3 && abs(cpu_result - 2.0f) < 1e-3;
    cout << "GPU梯形法积分结果：" << gpu_result << "，耗时：" << gpu_time << " ms" << endl;
    cout << "结果验证：" << (valid ? "正确" : "错误") << endl;
    cout << "GPU加速比：" << (double)cpu_time / gpu_time << "x" << endl;

    // 6. 释放资源
    CHECK_CUDA_ERR(cudaFree(d_samples));
    CHECK_CUDA_ERR(cudaFree(d_partial_sums));
    CHECK_CUDA_ERR(cudaEventDestroy(gpu_start));
    CHECK_CUDA_ERR(cudaEventDestroy(gpu_stop));

    return 0;
}

#include <cuda_runtime.h>
#include <iostream>
#include <vector>
#include <chrono>
#include <algorithm>

using namespace std;
using namespace chrono;

#define CHECK_CUDA_ERR(err) \
    if (err != cudaSuccess) { \
        cerr << "CUDA错误：" << cudaGetErrorString(err) << "（行号：" << __LINE__ << "）" << endl; \
        exit(1); \
    }

// 排序数组大小
const int ARRAY_SIZE = 100'000; // 10万个元素

// CPU串行排序（std::sort基准）
void sortCPU(vector<float>& arr) {
    sort(arr.begin(), arr.end());
}

// GPU核函数：奇偶排序单阶段（奇相位/偶相位）
__global__ void oddEvenSortStepGPU(float* d_arr, int size, bool* d_swapped, bool is_odd_phase) {
    int global_idx = threadIdx.x + blockIdx.x * blockDim.x;
    // 计算当前阶段的比较起始索引
    int i = is_odd_phase ? (2 * global_idx + 1) : (2 * global_idx);

    // 检查索引范围，比较并交换
    if (i < size - 1) {
        if (d_arr[i] > d_arr[i + 1]) {
            // 交换元素
            float temp = d_arr[i];
            d_arr[i] = d_arr[i + 1];
            d_arr[i + 1] = temp;
            *d_swapped = true; // 标记有交换
        }
    }
}

// GPU并行奇偶排序（主机端控制相位）
void oddEvenSortGPU(float* h_arr, int size) {
    float *d_arr;
    bool *d_swapped, h_swapped;
    const int BLOCK_SIZE = 256;
    const int GRID_SIZE = (size + BLOCK_SIZE - 1) / BLOCK_SIZE;

    // 设备内存分配
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_arr, size * sizeof(float)));
    CHECK_CUDA_ERR(cudaMalloc((void**)&d_swapped, sizeof(bool)));

    // 数据拷贝（主机→设备）
    CHECK_CUDA_ERR(cudaMemcpy(d_arr, h_arr, size * sizeof(float), cudaMemcpyHostToDevice));

    // 奇偶相位交替，直到无交换
    do {
        h_swapped = false;
        // 重置swapped标志（主机→设备）
        CHECK_CUDA_ERR(cudaMemcpy(d_swapped, &h_swapped, sizeof(bool), cudaMemcpyHostToDevice));

        // 1. 奇相位：比较(1,2), (3,4)...
        oddEvenSortStepGPU<<<GRID_SIZE, BLOCK_SIZE>>>(d_arr, size, d_swapped, true);
        CHECK_CUDA_ERR(cudaGetLastError());
        CHECK_CUDA_ERR(cudaDeviceSynchronize()); // 等待奇相位完成

        // 2. 偶相位：比较(0,1), (2,3)...
        oddEvenSortStepGPU<<<GRID_SIZE, BLOCK_SIZE>>>(d_arr, size, d_swapped, false);
        CHECK_CUDA_ERR(cudaGetLastError());
        CHECK_CUDA_ERR(cudaDeviceSynchronize()); // 等待偶相位完成

        // 读取swapped标志（设备→主机）
        CHECK_CUDA_ERR(cudaMemcpy(&h_swapped, d_swapped, sizeof(bool), cudaMemcpyDeviceToHost));
    } while (h_swapped);

    // 结果拷贝（设备→主机）
    CHECK_CUDA_ERR(cudaMemcpy(h_arr, d_arr, size * sizeof(float), cudaMemcpyDeviceToHost));

    // 释放资源
    CHECK_CUDA_ERR(cudaFree(d_arr));
    CHECK_CUDA_ERR(cudaFree(d_swapped));
}

int main() {
    // 1. 生成随机数组（CPU和GPU使用相同初始数据）
    vector<float> h_arr_cpu(ARRAY_SIZE);
    vector<float> h_arr_gpu(ARRAY_SIZE);
    srand(time(0));
    for (int i = 0; i < ARRAY_SIZE; ++i) {
        h_arr_cpu[i] = rand() / (float)RAND_MAX * 1000.0f; // 0~1000的随机数
        h_arr_gpu[i] = h_arr_cpu[i];
    }

    // 2. CPU串行排序
    auto cpu_start = high_resolution_clock::now();
    sortCPU(h_arr_cpu);
    auto cpu_time = duration_cast<milliseconds>(high_resolution_clock::now() - cpu_start).count();
    cout << "CPU std::sort耗时：" << cpu_time << " ms" << endl;

    // 3. GPU并行奇偶排序
    auto gpu_start = high_resolution_clock::now();
    oddEvenSortGPU(h_arr_gpu.data(), ARRAY_SIZE);
    auto gpu_time = duration_cast<milliseconds>(high_resolution_clock::now() - gpu_start).count();
    cout << "GPU奇偶排序耗时：" << gpu_time << " ms" << endl;

    // 4. 结果验证（对比前10个和后10个元素）
    bool valid = true;
    for (int i = 0; i < 10; ++i) {
        if (abs(h_arr_cpu[i] - h_arr_gpu[i]) > 1e-5) valid = false;
        if (abs(h_arr_cpu[ARRAY_SIZE - 10 + i] - h_arr_gpu[ARRAY_SIZE - 10 + i]) > 1e-5) valid = false;
    }
    cout << "结果验证：" << (valid ? "正确" : "错误") << endl;
    cout << "GPU加速比：" << (double)cpu_time / gpu_time << "x" << endl;

    return 0;
}