CUDAで配列の総和を求めたい

Question

前提

CUDAとCで非常に大きな一次元配列の総和を求めたいです。
double[100000000]のようなサイズです。

総和を求めるプログラムは下記URLをほぼそのまま使っています。
https://qiita.com/gyu-don/items/ef8a128fa24f6bddd342

下記コードだと配列数1048576までは総和を求められますが、
それ以上のサイズだと求めることができません。

大きな配列の総和を求めるようにするにはどのようにしたらよいのでしょうか。

CUDA(GPU)自体あまり理解できておらず、このような用途ではGPUはつかわないのか、
総和を求めるだけならCPUが早いのかなど一般的なCUDAの利用についても教えていただけるとたすかります。

該当のソースコード

#include <stdio.h>

__global__ void reduce0(double *g_idata, double *g_odata, unsigned int n){

        extern __shared__ double sdata[];

        unsigned int tid = threadIdx.x;
        unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;

        sdata[tid] = (i < n) ? g_idata[i] : 0;
        __syncthreads();

        for (unsigned int s=1; s<blockDim.x; s*=2) {
                if (tid % (2*s) == 0) {
                        sdata[tid] += sdata[tid + s];
                }
        __syncthreads();
        }

        if (tid == 0) g_odata[blockIdx.x] = sdata[0];

}

int main(int argc, char *argv[]){

        int threads_per_block = 1024;
        int number_of_blocks = 1024;

        //NG
        unsigned int n = 1048577;
        //OK
        //unsigned int n = 1048576;
        double array = n * sizeof(double);
        double result = number_of_blocks * sizeof(double);
        double *pt_array, *pt_result, cpu_result = 0.0, gpu_result = 0.0;
        cudaMallocManaged(&pt_array, array);
        cudaMallocManaged(&pt_result, result);

        // Initialize
        for(unsigned int i=0; i<n; i++){
                pt_array[i] = 2.0;
        }

        // CPU
        for(unsigned int i=0; i<n; i++){
                cpu_result = cpu_result + pt_array[i];
        }

        printf("cpu_result: %.1f\n", cpu_result);

        cudaFuncSetAttribute(&reduce0, cudaFuncAttributeMaxDynamicSharedMemorySize, 163840);
        reduce0<<< number_of_blocks, threads_per_block, 163840 >>>(pt_array, pt_result, n);
        cudaDeviceSynchronize();

        for(int i=0; i<number_of_blocks; i++){
                gpu_result = gpu_result + pt_result[i];
        }
        // GPU
        printf("gpu_result: %.1f\n", gpu_result);

        cudaFree(pt_array);
        cudaFree(pt_result);

}

試したこと

cudaFuncSetAttributeで共有メモリのサイズを拡張しました。

補足情報（FW/ツールのバージョンなど）

CUDA11.4で動作検証をしています。

Comments

[episteme]

> 総和を求めるだけならCPUが早いのか

CPU/GPUのふたつでやってみればいいじゃん

Answer 1

number_of_blocksを相当に増やせばいい、というか、nとthreads_per_block から計算するのが筋。

int number_of_blocks = ( n + threads_per_block  -1 ) / threads_per_block ;