GPU大设备数组优化：解决缓冲区不足及并行架构改进问询

问题背景 #

处理大规模设备数组时，设备内动态分配数组的dmemd内核在g_size≈8时出现设备缓冲区空间不足；改用主机侧分配数组的dmemh内核后，支持的g_size提升至约55k，但仍无法满足将g_size提升至120k的需求。

算法需求 #

• 完成g_size=120k规模的求和计算
• 在g_size循环内计算大型det数组
• 对det数组的特定索引进行求和
• 附示例代码为原求和逻辑的简化版本

当前代码架构 #

• 将串行代码中的g_size循环转为GPU上的g_size个线程块（CUDA最大块数为2^31）
• 每个线程块基于blockIdx创建det数组，再对det数组的特定索引求和
• 汇总所有线程块的结果得到最终总和

优化方案与并行化建议 #

1. 彻底消除det数组内存占用（最优方案） #

观察求和逻辑，仅需要det数组前g_sum_size个元素的和，完全可以在计算det元素的同时直接累加求和，无需存储整个det数组，彻底解决内存不足问题，同时减少内存读写操作提升效率。

优化后内核示例：

__global__ void no_det_kernel(cuDoubleComplex *d_sum_batch, int batch_start_idx) {
    int bid = blockIdx.x;
    int gid = batch_start_idx + bid;
    cuDoubleComplex dt_sum = make_cuDoubleComplex(0.0, 0.0);

    for (int i = 0; i < g_det_size; i++) {
        cuDoubleComplex val = make_cuDoubleComplex((i * 1.0 / (gid + 1)) * (0.01 / g_det_size),
                                      (i * 1.0 / (gid + 1)) * (0.01 / g_det_size));
        if (i < g_sum_size) {
            dt_sum = cuCadd(dt_sum, val);
        }
    }

    d_sum_batch[bid] = dt_sum;
}

2. 分批次处理线程块（兼容原有逻辑） #

若必须保留det数组的存储，可通过分批次处理线程块复用内存：每次仅处理N个线程块，复用同一块全局内存存储det数组，计算完一批后再处理下一批，最后汇总所有批次结果。

核心思路：

• 计算GPU能容纳的最大批次数量batch_size（需预留d_sum等其他内存空间）
• 循环处理每一批次，启动batch_size个线程块，处理完后将该批次的d_sum结果拷贝到主机临时存储，再复用d_det内存处理下一批
• 主机侧汇总所有批次的结果

3. 引入多线程并行提升计算效率 #

当前每个线程块仅使用1个线程，完全未发挥GPU多线程优势。可重构内核，让每个线程块包含多个线程，并行计算det数组元素，同时利用共享内存完成求和归约，减少全局内存访问延迟。

示例多线程内核：

__global__ void multi_thread_kernel(cuDoubleComplex *d_sum, cuDoubleComplex *d_det, int batch_start_idx) {
    int tid = threadIdx.x;
    int bid = blockIdx.x;
    int gid = batch_start_idx + bid;

    // 并行计算det数组元素，每个线程负责多个元素
    for (int i = tid; i < g_det_size; i += blockDim.x) {
        d_det[bid * g_det_size + i] = make_cuDoubleComplex((i * 1.0 / (gid + 1)) * (0.01 / g_det_size),
                                      (i * 1.0 / (gid + 1)) * (0.01 / g_det_size));
    }
    __syncthreads();

    // 并行计算局部和
    cuDoubleComplex local_sum = make_cuDoubleComplex(0.0, 0.0);
    for (int i = tid; i < g_sum_size; i += blockDim.x) {
        local_sum = cuCadd(local_sum, d_det[bid * g_det_size + i]);
    }

    // 共享内存归约求和
    __shared__ cuDoubleComplex s_sum[256];
    s_sum[tid] = local_sum;
    __syncthreads();

    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            s_sum[tid] = cuCadd(s_sum[tid], s_sum[tid + s]);
        }
        __syncthreads();
    }

    if (tid == 0) {
        d_sum[gid] = s_sum[0];
    }
}

完整优化后代码示例（无det数组版本） #

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <complex.h>
#include <cuComplex.h>

#include <math.h>
#include <string.h>
#include <time.h>
#include <ctime>
#include <stdint.h>
#include <cstring>

#include"cuda_runtime.h"
#include"device_launch_parameters.h"

// size of different arrays
const uint32_t g_target_size = 120000U;
const uint32_t g_det_size = 40000U;
const uint32_t g_sum_size = 30000U;

// 优化内核：无det数组直接求和
__global__ void no_det_kernel(cuDoubleComplex *d_sum_batch, int batch_start_idx) {
    int bid = blockIdx.x;
    int gid = batch_start_idx + bid;
    cuDoubleComplex dt_sum = make_cuDoubleComplex(0.0, 0.0);

    for (int i = 0; i < g_det_size; i++) {
        cuDoubleComplex val = make_cuDoubleComplex((i * 1.0 / (gid + 1)) * (0.01 / g_det_size),
                                      (i * 1.0 / (gid + 1)) * (0.01 / g_det_size));
        if (i < g_sum_size) {
            dt_sum = cuCadd(dt_sum, val);
        }
    }

    d_sum_batch[bid] = dt_sum;
}

int main() {
    double complex tsum = 0.0 + 0.0 * I;

    // 分批次处理，可根据GPU显存调整batch_size
    const int batch_size = 10240;
    int total_batches = (g_target_size + batch_size - 1) / batch_size;

    // 分配设备侧批次求和结果内存
    cuDoubleComplex *d_sum_batch;
    cudaMalloc((void **)&d_sum_batch, sizeof(cuDoubleComplex) * batch_size);
    // 主机侧临时存储批次结果
    double complex *sum_batch = (double complex*)malloc(sizeof(double complex) * batch_size);

    for (int batch_idx = 0; batch_idx < total_batches; batch_idx++) {
        int current_batch_size = (batch_idx == total_batches - 1) ? (g_target_size - batch_idx * batch_size) : batch_size;
        int start_idx = batch_idx * batch_size;

        dim3 block(1); // 可改为block(256)进一步优化计算并行度
        dim3 grid(current_batch_size);

        no_det_kernel <<< grid, block >>> (d_sum_batch, start_idx);
        cudaDeviceSynchronize();

        // 拷贝批次结果到主机
        cudaMemcpy(sum_batch, d_sum_batch, sizeof(double complex) * current_batch_size, cudaMemcpyDeviceToHost);

        // 累加批次结果
        for (int i = 0; i < current_batch_size; i++) {
            tsum += sum_batch[i];
        }
    }

    printf("Total sum: %.16E,%.16E \n", creal(tsum), cimag(tsum));

    // 释放内存
    free(sum_batch);
    cudaFree(d_sum_batch);
    cudaDeviceReset();

    return 0;
}

内容的提问来源于stack exchange，提问作者Anomalous Physicist